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

Merge remote-tracking branch 'upstream/develop' into outer-product

Browse source
This commit is contained in:
Sam Shaner 2015-11-25 13:09:11 -08:00
commit 74ffcb4475
27 changed files with 650 additions and 938 deletions
Download patch file Download diff file Expand all files Collapse all files

2
openmc/mgxs/library.py
View file

@ -361,6 +361,8 @@ class Library(object):
raise ValueError(msg)
self._sp_filename = statepoint._f.filename
self._openmc_geometry = statepoint.summary.openmc_geometry
self._opencg_geometry = None
# Load tallies for each MGXS for each domain and mgxs type
for domain in self.domains:

2
openmc/opencg_compatible.py
View file

@ -726,7 +726,7 @@ def get_openmc_cell(opencg_cell):
openmc_cell.fill = get_openmc_material(fill)
if opencg_cell.rotation:
rotation = np.asarray(opencg_cell.rotation, dtype=np.int)
rotation = np.asarray(opencg_cell.rotation, dtype=np.float64)
openmc_cell.rotation = rotation
if opencg_cell.translation:

76
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
@ -122,7 +124,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
@ -196,9 +200,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
@ -399,6 +407,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)
@ -861,57 +877,47 @@ 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)
@ -1182,10 +1188,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_source_subelement()
self._create_output_subelement()
self._create_statepoint_subelement()

404
src/ace.F90
View file

@ -45,9 +45,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), pointer :: nuc => null()
type(SAlphaBeta), pointer :: sab => null()
type(Material), pointer :: mat
type(Nuclide), pointer :: nuc
type(SAlphaBeta), pointer :: sab
type(SetChar) :: already_read
! allocate arrays for ACE table storage and cross section cache
@ -234,7 +234,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
@ -258,9 +257,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), pointer :: nuc => null()
type(SAlphaBeta), pointer :: sab => null()
type(XsListing), pointer :: listing => null()
type(Nuclide), pointer :: nuc
type(SAlphaBeta), pointer :: sab
type(XsListing), pointer :: listing
! determine path, record length, and location of table
listing => xs_listings(i_listing)
@ -406,8 +405,6 @@ contains
end select
deallocate(XSS)
if(associated(nuc)) nullify(nuc)
if(associated(sab)) nullify(sab)
end subroutine read_ace_table
@ -417,10 +414,8 @@ contains
!===============================================================================
subroutine read_esz(nuc, data_0K)
type(Nuclide), pointer :: nuc
logical :: data_0K ! are we reading 0K data?
type(Nuclide), 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
@ -507,8 +502,7 @@ contains
!===============================================================================
subroutine read_nu_data(nuc)
type(Nuclide), pointer :: nuc
type(Nuclide), intent(inout) :: nuc
integer :: i ! loop index
integer :: JXS2 ! location for fission nu data
@ -524,7 +518,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)
@ -707,8 +701,7 @@ contains
!===============================================================================
subroutine read_reactions(nuc)
type(Nuclide), pointer :: nuc
type(Nuclide), intent(inout) :: nuc
integer :: i ! loop indices
integer :: i_fission ! index in nuc % index_fission
@ -722,7 +715,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)
@ -740,14 +732,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
@ -760,123 +753,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
@ -890,8 +885,7 @@ contains
!===============================================================================
subroutine read_angular_dist(nuc)
type(Nuclide), pointer :: nuc
type(Nuclide), intent(inout) :: nuc
integer :: JXS8 ! location of angular distribution locators
integer :: JXS9 ! location of angular distributions
@ -902,7 +896,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)
@ -910,71 +903,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
@ -985,29 +979,28 @@ contains
!===============================================================================
subroutine read_energy_dist(nuc)
type(Nuclide), pointer :: nuc
type(Nuclide), 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
@ -1019,10 +1012,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
@ -1102,7 +1094,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
@ -1146,7 +1137,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
@ -1204,7 +1195,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
@ -1234,7 +1225,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
@ -1285,7 +1276,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
@ -1301,8 +1292,7 @@ contains
!===============================================================================
subroutine read_unr_res(nuc)
type(Nuclide), pointer :: nuc
type(Nuclide), intent(inout) :: nuc
integer :: JXS23 ! location of URR data
integer :: lc ! locator
@ -1390,8 +1380,7 @@ contains
!===============================================================================
subroutine generate_nu_fission(nuc)
type(Nuclide), pointer :: nuc
type(Nuclide), intent(inout) :: nuc
integer :: i ! index on nuclide energy grid
real(8) :: E ! energy
@ -1417,8 +1406,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
@ -1600,7 +1588,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

96
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,10 +76,6 @@ 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
!===============================================================================
@ -99,7 +92,7 @@ module ace_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
! Energy grid information
integer :: n_grid ! # of nuclide grid points
@ -153,7 +146,9 @@ module ace_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
@ -166,14 +161,12 @@ module ace_header
!===============================================================================
type Nuclide0K
character(10) :: nuclide ! name of nuclide, e.g. U-238
character(16) :: scheme = 'ares' ! target velocity sampling scheme
character(10) :: name ! name of nuclide, e.g. 92235.03c
character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c
real(8) :: E_min = 0.01e-6_8 ! lower cutoff energy for res scattering
real(8) :: E_max = 1000.0e-6_8 ! upper cutoff energy for res scattering
end type Nuclide0K
!===============================================================================
@ -265,7 +258,6 @@ module ace_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)
@ -288,24 +280,10 @@ module ace_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
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.
!===============================================================================
@ -314,12 +292,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()
@ -336,32 +308,13 @@ 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
!===============================================================================
! NUCLIDE_CLEAR resets and deallocates data in Nuclide.
!===============================================================================
@ -372,31 +325,6 @@ module ace_header
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()
@ -405,18 +333,16 @@ module ace_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 this % nuc_list % clear()
call this % reaction_index % clear()
end subroutine nuclide_clear

95
src/cross_section.F90
View file

@ -13,10 +13,6 @@ module cross_section
use search, only: binary_search
implicit none
save
integer :: union_grid_index
!$omp threadprivate(union_grid_index)
contains
@ -33,7 +29,8 @@ contains
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
@ -44,7 +41,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
@ -52,11 +48,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
@ -99,9 +94,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
! ========================================================================
@ -129,10 +124,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
@ -142,18 +133,19 @@ 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
real(8) :: f ! interp factor on nuclide energy grid
type(Nuclide), pointer :: nuc
type(Material), pointer :: mat
@ -165,7 +157,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
@ -178,8 +170,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), &
@ -219,7 +211,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) &
@ -241,13 +232,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,7 +366,6 @@ contains
logical :: same_nuc ! do we know the xs for this nuclide at this energy?
type(UrrData), pointer :: urr
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
micro_xs(i_nuclide) % use_ptable = .true.
@ -408,7 +391,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
@ -416,7 +399,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
@ -478,18 +461,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
@ -526,38 +508,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), intent(in) :: nuc ! target nuclide at temperature
real(8) :: xs_out ! 0K xs at trial energy
type(Nuclide), 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

28
src/fission.F90
View file

@ -15,18 +15,17 @@ contains
! given nuclide and incoming neutron energy
!===============================================================================
function nu_total(nuc, E) result(nu)
type(Nuclide), pointer :: 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
pure function nu_total(nuc, E) result(nu)
type(Nuclide), 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
integer :: NC ! number of polynomial coefficients
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,11 +48,10 @@ contains
! for a given nuclide and incoming neutron energy
!===============================================================================
function nu_prompt(nuc, E) result(nu)
type(Nuclide), pointer :: 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
pure function nu_prompt(nuc, E) result(nu)
type(Nuclide), 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
integer :: NC ! number of polynomial coefficients
@ -87,10 +85,9 @@ contains
! for a given nuclide and incoming neutron energy
!===============================================================================
function nu_delayed(nuc, E) result(nu)
pure function nu_delayed(nuc, E) result(nu)
type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
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
@ -111,8 +108,7 @@ contains
! a given nuclide and incoming neutron energy in a given delayed group.
!===============================================================================
function yield_delayed(nuc, E, g) result(yield)
pure function yield_delayed(nuc, E, g) result(yield)
type(Nuclide), 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

89
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 string, only: to_str
use tally, only: score_surface_current
@ -871,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
!===============================================================================

16
src/global.F90
View file

@ -436,7 +436,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
@ -463,14 +468,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)

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/math.F90
View file

@ -12,7 +12,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
@ -71,7 +71,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
@ -123,7 +123,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

@ -18,9 +18,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
@ -71,9 +70,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?
@ -96,11 +94,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
@ -130,9 +127,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(:)
@ -167,9 +163,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
@ -262,9 +258,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
@ -328,9 +323,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

99
src/output.F90
View file

@ -100,10 +100,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
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
@ -195,9 +194,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
@ -250,7 +248,6 @@ contains
!===============================================================================
subroutine print_particle(p)
type(Particle), intent(in) :: p
integer :: i ! index for coordinate levels
@ -320,9 +317,8 @@ contains
!===============================================================================
subroutine print_nuclide(nuc, unit)
type(Nuclide), pointer :: nuc
integer, optional :: unit
type(Nuclide), intent(in) :: nuc
integer, intent(in), optional :: unit
integer :: i ! loop index over nuclides
integer :: unit_ ! unit to write to
@ -334,8 +330,7 @@ contains
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
@ -363,32 +358,32 @@ contains
! Information on each reaction
write(unit_,*) ' Reaction Q-value COM Law IE size(angle) size(energy)'
do i = 1, nuc % n_reaction
rxn => nuc % reactions(i)
associate (rxn => nuc % 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 + (nuc % 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 + (nuc % 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,
@ -438,9 +433,8 @@ contains
!===============================================================================
subroutine print_sab_table(sab, unit)
type(SAlphaBeta), pointer :: sab
integer, optional :: unit
type(SAlphaBeta), intent(in) :: sab
integer, intent(in), optional :: unit
integer :: size_sab ! memory used by S(a,b) table
integer :: unit_ ! unit to write to
@ -526,8 +520,8 @@ contains
integer :: i ! loop index
integer :: unit_xs ! cross_sections.out file unit
character(MAX_FILE_LEN) :: path ! path of summary file
type(Nuclide), pointer :: nuc => null()
type(SAlphaBeta), pointer :: sab => null()
type(Nuclide), pointer :: nuc
type(SAlphaBeta), pointer :: sab
! Create filename for log file
path = trim(path_output) // "cross_sections.out"
@ -681,7 +675,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")
@ -1183,7 +1177,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
@ -1355,10 +1349,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
@ -1422,12 +1415,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
@ -1439,7 +1432,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

86
src/physics.F90
View file

@ -185,7 +185,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
@ -195,7 +194,6 @@ contains
real(8) :: prob
real(8) :: cutoff
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
! Get pointer to nuclide
nuc => nuclides(i_nuclide)
@ -220,14 +218,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
@ -312,11 +311,10 @@ contains
real(8) :: f
real(8) :: prob
real(8) :: cutoff
type(Nuclide), 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), pointer :: nuc
! copy incoming direction
uvw_old(:) = p % coord(1) % uvw
@ -343,11 +341,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
@ -370,28 +365,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
@ -420,9 +415,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
@ -759,9 +753,7 @@ contains
!===============================================================================
subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff)
type(Nuclide), pointer :: nuc ! target nuclide at temperature T
type(Nuclide), 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
@ -1006,8 +998,7 @@ contains
!===============================================================================
subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw)
type(Nuclide), pointer :: nuc ! target nuclide at temperature
type(Nuclide), 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)
@ -1080,7 +1071,6 @@ 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
@ -1095,11 +1085,9 @@ contains
real(8) :: weight ! weight adjustment for ufs method
logical :: in_mesh ! source site in ufs mesh?
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
! Get pointers
nuc => nuclides(i_nuclide)
rxn => nuc % reactions(i_reaction)
! TODO: Heat generation from fission
@ -1170,7 +1158,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
@ -1197,8 +1186,8 @@ contains
function sample_fission_energy(nuc, rxn, p) result(E_out)
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
type(Nuclide), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
type(Particle), intent(inout) :: p ! Particle causing fission
real(8) :: E_out ! outgoing energy of fission neutron
@ -1323,8 +1312,8 @@ contains
!===============================================================================
subroutine inelastic_scatter(nuc, rxn, p)
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
type(Nuclide), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
type(Particle), intent(inout) :: p
integer :: i ! loop index
@ -1409,8 +1398,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)
@ -1536,7 +1524,6 @@ contains
!===============================================================================
function rotate_angle(uvw0, mu) result(uvw)
real(8), intent(in) :: uvw0(3) ! directional cosine
real(8), intent(in) :: mu ! cosine of angle in lab or CM
real(8) :: uvw(3) ! rotated directional cosine
@ -1585,8 +1572,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

3
src/source.F90
View file

@ -96,8 +96,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

82
src/string.F90
View file

@ -25,7 +25,6 @@ contains
!===============================================================================
subroutine split_string(string, words, n)
character(*), intent(in) :: string
character(*), intent(out) :: words(MAX_WORDS)
integer, intent(out) :: n
@ -166,7 +165,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)
@ -186,8 +185,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
@ -209,8 +207,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
@ -234,39 +231,38 @@ contains
! integers.
!===============================================================================
function zero_padded(num, n_digits) result(str)
integer, intent(in) :: num
integer, intent(in) :: n_digits
character(11) :: str
function zero_padded(num, n_digits) result(str)
integer, intent(in) :: num
integer, intent(in) :: n_digits
character(11) :: str
character(8) :: zp_form
character(8) :: zp_form
! Make sure n_digits is reasonable. 10 digits is the maximum needed for the
! largest integer(4).
if (n_digits > 10) then
call fatal_error('zero_padded called with an unreasonably large &
&n_digits (>10)')
end if
! Make sure n_digits is reasonable. 10 digits is the maximum needed for the
! largest integer(4).
if (n_digits > 10) then
call fatal_error('zero_padded called with an unreasonably large &
&n_digits (>10)')
end if
! Write a format string of the form '(In.m)' where n is the max width and
! m is the min width. If a sign is present, then n must be one greater
! than m.
if (num < 0) then
write(zp_form, '("(I", I0, ".", I0, ")")') n_digits+1, n_digits
else
write(zp_form, '("(I", I0, ".", I0, ")")') n_digits, n_digits
end if
! Write a format string of the form '(In.m)' where n is the max width and
! m is the min width. If a sign is present, then n must be one greater
! than m.
if (num < 0) then
write(zp_form, '("(I", I0, ".", I0, ")")') n_digits+1, n_digits
else
write(zp_form, '("(I", I0, ".", I0, ")")') n_digits, n_digits
end if
! Format the number.
write(str, zp_form) num
end function zero_padded
! Format the number.
write(str, zp_form) num
end function zero_padded
!===============================================================================
! 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
@ -286,10 +282,9 @@ end function zero_padded
! 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
@ -321,10 +316,9 @@ end function zero_padded
! 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
@ -350,7 +344,7 @@ end function zero_padded
! integer.
!===============================================================================
function count_digits(num) result(n_digits)
pure function count_digits(num) result(n_digits)
integer, intent(in) :: num
integer :: n_digits
@ -368,7 +362,7 @@ end function zero_padded
! 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
@ -382,7 +376,7 @@ end function zero_padded
! 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
@ -396,7 +390,7 @@ end function zero_padded
! 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
@ -421,7 +415,7 @@ end function zero_padded
! 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
@ -440,7 +434,7 @@ end function zero_padded
! 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

252
src/tally.F90
View file

@ -37,13 +37,13 @@ contains
subroutine score_general(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
@ -66,9 +66,6 @@ contains
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), pointer :: nuc
i = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
@ -220,24 +217,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
@ -257,24 +250,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
@ -294,24 +283,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
@ -466,9 +451,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
@ -480,11 +462,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
@ -494,7 +476,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
@ -535,9 +517,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
@ -548,11 +527,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
@ -560,27 +539,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
@ -588,15 +564,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
@ -606,13 +579,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 &
@ -624,38 +597,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
@ -691,14 +693,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
@ -708,26 +706,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
@ -737,9 +729,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
@ -1018,9 +1009,8 @@ contains
!===============================================================================
subroutine score_fission_eout(p, t, i_score)
type(Particle), intent(in) :: p
type(TallyObject), pointer :: t
type(TallyObject), intent(inout) :: t
integer, intent(in) :: i_score ! index for score
integer :: i ! index of outgoing energy filter
@ -1348,9 +1338,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
@ -1724,7 +1714,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.
@ -1948,7 +1938,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
!===============================================================================
@ -129,81 +125,6 @@ module tally_header
! Tally precision triggers
integer :: n_triggers = 0 ! # of triggers
type(TriggerObject), allocatable :: triggers(:) ! Array of triggers
! 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