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

Merge remote-tracking branch 'samuelshaner/photon-new' into photon-new

Browse source
This commit is contained in:
Paul Romano 2017-07-19 11:08:26 -05:00
commit cd185bb33d
16 changed files with 755 additions and 207 deletions
Download patch file Download diff file Expand all files Collapse all files

23
openmc/data/library.py
View file

@ -63,32 +63,29 @@ class DataLibrary(EqualityMixin):
library = {'path': filename, 'type': filetype, 'materials': materials}
self.libraries.append(library)
def export_to_xml(self, path='cross_sections.xml'):
def export_to_xml(self, path='cross_sections.xml', append=False):
"""Export cross section data library to an XML file.
Parameters
----------
path : str
Path to file to write. Defaults to 'cross_sections.xml'.
append : bool
Whether to append to an existing file, if it exists.
Defaults to False.
"""
root = ET.Element('cross_sections')
# Determine common directory for library paths
common_dir = os.path.dirname(os.path.commonprefix(
[lib['path'] for lib in self.libraries]))
if common_dir == '':
common_dir = '.'
directory = os.path.relpath(common_dir, os.path.dirname(path))
if directory != '.':
dir_element = ET.SubElement(root, "directory")
dir_element.text = directory
if append:
root = ET.parse(path).getroot()
else:
root = ET.Element('cross_sections')
for library in self.libraries:
lib_element = ET.SubElement(root, "library")
lib_element.set('materials', ' '.join(library['materials']))
lib_element.set('path', os.path.relpath(library['path'], common_dir))
lib_element.set('path', os.path.relpath(library['path'],
os.path.dirname(path)))
lib_element.set('type', library['type'])
# Clean the indentation to be user-readable

40
openmc/filter.py
View file

@ -17,7 +17,7 @@ from .mixin import IDManagerMixin
_FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface',
'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal',
'distribcell', 'delayedgroup', 'energyfunction']
'distribcell', 'delayedgroup', 'energyfunction', 'particle']
_CURRENT_NAMES = {1: 'x-min out', 2: 'x-min in',
3: 'x-max out', 4: 'x-max in',
@ -26,6 +26,7 @@ _CURRENT_NAMES = {1: 'x-min out', 2: 'x-min in',
9: 'z-min out', 10: 'z-min in',
11: 'z-max out', 12: 'z-max in'}
_PARTICLE_IDS = {'neutron': 1, 'photon': 2, 'electron': 3, 'positron': 4}
class FilterMeta(ABCMeta):
def __new__(cls, name, bases, namespace, **kwargs):
@ -549,6 +550,43 @@ class MaterialFilter(WithIDFilter):
self._smart_set_bins(bins, openmc.Material)
class ParticleFilter(WithIDFilter):
"""Bins tally event locations based on the Particle type.
Parameters
----------
bins : Str, Integral, or iterable thereof
The Particles to tally. Either str with particle type or their
Integral ID numbers can be used with IDs listed in _PARTICLE_IDS.
filter_id : int
Unique identifier for the filter
Attributes
----------
bins : Iterable of Integral
openmc.Materi IDs.
id : int
Unique identifier for the filter
num_bins : Integral
The number of filter bins
stride : Integral
The number of filter, nuclide and score bins within each of this
filter's bins.
"""
@property
def bins(self):
return self._bins
@bins.setter
def bins(self, bins):
bins = np.atleast_1d(bins)
cv.check_iterable_type('filter bins', bins, str)
bins = np.atleast_1d([b if isinstance(b, Integral) else _PARTICLE_IDS[b]
for b in bins])
self._bins = bins
class CellFilter(WithIDFilter):
"""Bins tally event locations based on the Cell they occured in.

26
openmc/settings.py
View file

@ -38,12 +38,13 @@ class Settings(object):
calculations to find the path to the XML cross section file.
cutoff : dict
Dictionary defining weight cutoff and energy cutoff. The dictionary may
have three keys, 'weight', 'weight_avg' and 'energy'. Value for 'weight'
have six keys, 'weight', 'weight_avg', 'energy_neutron', 'energy_photon',
'energy_electron', and 'energy_positron'. Value for 'weight'
should be a float indicating weight cutoff below which particle undergo
Russian roulette. Value for 'weight_avg' should be a float indicating
weight assigned to particles that are not killed after Russian
roulette. Value of energy should be a float indicating energy in eV
below which particle will be killed.
below which particle type will be killed.
energy_mode : {'continuous-energy', 'multi-group'}
Set whether the calculation should be continuous-energy or multi-group.
entropy_mesh : openmc.Mesh
@ -79,6 +80,8 @@ class Settings(object):
:tallies: Whether the 'tallies.out' file should be written (bool)
particles : int
Number of particles per generation
photon_transport : bool
Whether to use photon transport.
ptables : bool
Determine whether probability tables are used.
resonance_scattering : dict
@ -183,6 +186,7 @@ class Settings(object):
self._confidence_intervals = None
self._cross_sections = None
self._multipole_library = None
self._photon_transport = None
self._ptables = None
self._run_cmfd = None
self._seed = None
@ -286,6 +290,10 @@ class Settings(object):
def ptables(self):
return self._ptables
@property
def photon_transport(self):
return self._photon_transport
@property
def run_cmfd(self):
return self._run_cmfd
@ -548,6 +556,11 @@ class Settings(object):
cv.check_type('multipole library', multipole_library, string_types)
self._multipole_library = multipole_library
@photon_transport.setter
def photon_transport(self, photon_transport):
cv.check_type('photon transport', photon_transport, bool)
self._photon_transport = photon_transport
@ptables.setter
def ptables(self, ptables):
cv.check_type('probability tables', ptables, bool)
@ -583,7 +596,8 @@ class Settings(object):
cv.check_type('average survival weight', cutoff[key], Real)
cv.check_greater_than('average survival weight',
cutoff[key], 0.0)
elif key in ['energy', 'energy_photon']:
elif key in ['energy_neutron', 'energy_photon', 'energy_electron',
'energy_positron']:
cv.check_type('energy cutoff', cutoff[key], Real)
cv.check_greater_than('energy cutoff', cutoff[key], 0.0)
else:
@ -922,6 +936,11 @@ class Settings(object):
element = ET.SubElement(root, "multipole_library")
element.text = str(self._multipole_library)
def _create_photon_transport_subelement(self, root):
if self._photon_transport is not None:
element = ET.SubElement(root, "photon_transport")
element.text = str(self._photon_transport).lower()
def _create_ptables_subelement(self, root):
if self._ptables is not None:
element = ET.SubElement(root, "ptables")
@ -1112,6 +1131,7 @@ class Settings(object):
self._create_multipole_library_subelement(root_element)
self._create_energy_mode_subelement(root_element)
self._create_max_order_subelement(root_element)
self._create_photon_transport_subelement(root_element)
self._create_ptables_subelement(root_element)
self._create_run_cmfd_subelement(root_element)
self._create_seed_subelement(root_element)

11
scripts/openmc-get-nndc-data
View file

@ -33,6 +33,8 @@ parser = argparse.ArgumentParser(
)
parser.add_argument('-b', '--batch', action='store_true',
help='supresses standard in')
parser.add_argument('-p', '--photo', default='generate_true',
help='Whether to include photo-atomic interaction data')
args = parser.parse_args()
@ -111,7 +113,7 @@ for f in files:
# Move ACE files down one level
for filename in glob.glob('nndc/293.6K/ENDF-B-VII.1-neutron-293.6K/*'):
shutil.move(filename, 'nndc/293.6K/')
shutil.move(filename, 'nndc/293.6K/' + os.path.basename(filename))
# ==============================================================================
# FIX ZAID ASSIGNMENTS FOR VARIOUS S(A,B) TABLES
@ -158,3 +160,10 @@ pwd = os.path.dirname(os.path.realpath(__file__))
ace2hdf5 = os.path.join(pwd, 'openmc-ace-to-hdf5')
subprocess.call([ace2hdf5, '-d', 'nndc_hdf5', '--fission_energy_release',
fer_file] + ace_files)
# Generate photo interaction library files
if args.photo == 'generate_true':
pwd = os.path.dirname(os.path.realpath(__file__))
photo_endf = os.path.join(pwd, 'openmc-get-photo-endf71')
subprocess.call([photo_endf, '-c', 'nndc_hdf5/cross_sections.xml'])

85
scripts/openmc-get-photo-endf71 Executable file
View file

@ -0,0 +1,85 @@
#!/usr/bin/env python
from __future__ import print_function
import os
import sys
import shutil
import zipfile
import requests
import argparse
from io import BytesIO
import openmc.data
from openmc.data import ATOMIC_SYMBOL
description = """
Download ENDF/B-VII.1 ENDF data from NNDC for photo-atomic and atomic
relaxation data and convert it to an HDF5 library for use with OpenMC.
This data is used for photon transport in OpenMC.
"""
class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter,
argparse.RawDescriptionHelpFormatter):
pass
parser = argparse.ArgumentParser(
description=description,
formatter_class=CustomFormatter
)
parser.add_argument('-c', '--cross-sections-file',
help='cross_sections.xml file to append libraries to')
args = parser.parse_args()
base_url = 'http://www.nndc.bnl.gov/endf/b7.1/zips/'
files = ['ENDF-B-VII.1-photoat.zip', 'ENDF-B-VII.1-atomic_relax.zip']
# ==============================================================================
# DOWNLOAD FILES FROM NNDC SITE
if not os.path.exists('photo_hdf5'):
os.mkdir('photo_hdf5')
library = openmc.data.DataLibrary()
filesComplete = []
for f in files:
# Establish connection to URL
print('Downloading {}...'.format(f))
url = base_url + f
r = requests.get(url, stream=True)
zipfile.ZipFile(BytesIO(r.content)).extractall()
# ==============================================================================
# GENERATE HDF5 LIBRARY
for z in range(1, 101):
element = ATOMIC_SYMBOL[z]
print('Extracting {} interaction data...'.format(element))
# Load files
filename = 'photoat/photoat-{:03}_{}_000.endf'.format(z, element)
photo_file = 'photoat/' + element + '.endf'
shutil.move(filename, photo_file)
filename = 'atomic_relax/atom-{:03}_{}_000.endf'.format(z, element)
atom_file = 'atomic_relax/' + element + '.endf'
shutil.move(filename, atom_file)
hdf5_file = 'photo_hdf5/' + element + '.h5'
if os.path.isfile(hdf5_file):
os.remove(hdf5_file)
f = openmc.data.IncidentPhoton.from_endf(photo_file, atom_file)
f.export_to_hdf5(hdf5_file)
library.register_file(hdf5_file)
if args.cross_sections_file is not None:
path = args.cross_sections_file
library.export_to_xml(path, True)
else:
path = 'photo_hdf5/cross_sections.xml'
library.export_to_xml(path)

20
src/constants.F90
View file

@ -179,7 +179,8 @@ module constants
integer, parameter :: &
NEUTRON = 1, &
PHOTON = 2, &
ELECTRON = 3
ELECTRON = 3, &
POSITRON = 4
! Angular distribution type
integer, parameter :: &
@ -340,8 +341,8 @@ module constants
SCORE_DELAYED_NU_FISSION = -19, & ! delayed neutron production rate
SCORE_PROMPT_NU_FISSION = -20, & ! prompt neutron production rate
SCORE_INVERSE_VELOCITY = -21, & ! flux-weighted inverse velocity
SCORE_FISS_Q_PROMPT = -22, & ! prompt fission Q-value
SCORE_FISS_Q_RECOV = -23, & ! recoverable fission Q-value
SCORE_FISS_Q_RECOV = -22, & ! recoverable fission Q-value
SCORE_FISS_Q_PROMPT = -23, & ! prompt fission Q-value
SCORE_DECAY_RATE = -24 ! delayed neutron precursor decay rate
! Maximum scattering order supported
@ -365,7 +366,7 @@ module constants
integer, parameter :: NO_BIN_FOUND = -1
! Tally filter and map types
integer, parameter :: N_FILTER_TYPES = 14
integer, parameter :: N_FILTER_TYPES = 15
integer, parameter :: &
FILTER_UNIVERSE = 1, &
FILTER_MATERIAL = 2, &
@ -380,7 +381,8 @@ module constants
FILTER_POLAR = 11, &
FILTER_AZIMUTHAL = 12, &
FILTER_DELAYEDGROUP = 13, &
FILTER_ENERGYFUNCTION = 14
FILTER_ENERGYFUNCTION = 14, &
FILTER_PARTICLE = 15
! Mesh types
integer, parameter :: &
@ -424,12 +426,13 @@ module constants
! ============================================================================
! RANDOM NUMBER STREAM CONSTANTS
integer, parameter :: N_STREAMS = 5
integer, parameter :: N_STREAMS = 6
integer, parameter :: STREAM_TRACKING = 1
integer, parameter :: STREAM_TALLIES = 2
integer, parameter :: STREAM_SOURCE = 3
integer, parameter :: STREAM_URR_PTABLE = 4
integer, parameter :: STREAM_VOLUME = 5
integer, parameter :: STREAM_PHOTON = 6
! ============================================================================
! MISCELLANEOUS CONSTANTS
@ -450,6 +453,11 @@ module constants
MODE_PARTICLE = 4, & ! Particle restart mode
MODE_VOLUME = 5 ! Volume calculation mode
! Electron treatments
integer, parameter :: &
ELECTRON_LED = 1, & ! Local Energy Deposition
ELECTRON_TTB = 2 ! Thick Target Bremsstrahlung
!=============================================================================
! CMFD CONSTANTS

35
src/cross_section.F90
View file

@ -27,6 +27,17 @@ contains
subroutine calculate_xs(p)
type(Particle), intent(inout) :: p
! Set all material macroscopic cross sections to zero
material_xs % total = ZERO
material_xs % elastic = ZERO
material_xs % absorption = ZERO
material_xs % fission = ZERO
material_xs % nu_fission = ZERO
material_xs % coherent = ZERO
material_xs % incoherent = ZERO
material_xs % photoelectric = ZERO
material_xs % pair_production = ZERO
if (p % type == NEUTRON) then
call calculate_neutron_xs(p)
elseif (p % type == PHOTON) then
@ -52,13 +63,6 @@ contains
real(8) :: atom_density ! atom density of a nuclide
logical :: check_sab ! should we check for S(a,b) table?
! Set all material macroscopic cross sections to zero
material_xs % total = ZERO
material_xs % elastic = ZERO
material_xs % absorption = ZERO
material_xs % fission = ZERO
material_xs % nu_fission = ZERO
! Exit subroutine if material is void
if (p % material == MATERIAL_VOID) return
@ -257,8 +261,9 @@ contains
micro_xs(i_nuclide) % interp_factor = f
! Initialize nuclide cross-sections to zero
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
micro_xs(i_nuclide) % photon_prod = ZERO
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % total = (ONE - f) * xs % total(i_grid) &
@ -272,6 +277,10 @@ contains
micro_xs(i_nuclide) % absorption = (ONE - f) * xs % absorption( &
i_grid) + f * xs % absorption(i_grid + 1)
! Calculate microscopic nuclide photon production cross section
micro_xs(i_nuclide) % photon_prod = (ONE - f) * xs % &
photon_prod(i_grid) + f * xs % photon_prod(i_grid + 1)
if (nuc % fissionable) then
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % fission = (ONE - f) * xs % fission(i_grid) &
@ -442,6 +451,7 @@ contains
integer :: i_energy ! index for energy
integer :: i_low ! band index at lower bounding energy
integer :: i_up ! band index at upper bounding energy
integer :: i_grid ! index on nuclide energy grid
real(8) :: f ! interpolation factor
real(8) :: r ! pseudo-random number
real(8) :: elastic ! elastic cross section
@ -584,13 +594,6 @@ contains
integer :: i_element ! index into elements array
real(8) :: atom_density ! atom density of a nuclide
! Set all material macroscopic cross sections to zero
material_xs % total = ZERO
material_xs % coherent = ZERO
material_xs % incoherent = ZERO
material_xs % photoelectric = ZERO
material_xs % pair_production = ZERO
! Exit subroutine if material is void
if (p % material == MATERIAL_VOID) return

6
src/energy_distribution.F90
View file

@ -333,8 +333,8 @@ contains
c_k = c_k1
end do
! Check to make sure k is <= NP - 1
k = min(k, n_energy_out - 1)
! Check to make sure 1 <= k <= NP - 1
k = max(1, min(k, n_energy_out - 1))
E_l_k = this%distribution(l)%e_out(k)
p_l_k = this%distribution(l)%p(k)
@ -361,7 +361,7 @@ contains
end if
! Now interpolate between incident energy bins i and i + 1
if (.not. histogram_interp) then
if (.not. histogram_interp .and. n_energy_out > 1) then
if (l == i) then
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
else

3
src/global.F90
View file

@ -121,6 +121,7 @@ module global
real(8) :: log_spacing ! spacing on logarithmic grid
logical :: photon_transport = .false.
integer :: electron_treatment = ELECTRON_LED
! ============================================================================
! MULTI-GROUP CROSS SECTION RELATED VARIABLES
@ -313,7 +314,7 @@ module global
logical :: survival_biasing = .false.
real(8) :: weight_cutoff = 0.25_8
real(8) :: energy_cutoff(3) = [ZERO, 1000.0_8, ZERO]
real(8) :: energy_cutoff(4) = [ZERO, 1000.0_8, ZERO, ZERO]
real(8) :: weight_survive = ONE
! ============================================================================

100
src/input_xml.F90
View file

@ -267,6 +267,30 @@ contains
! Copy random number seed if specified
if (check_for_node(root, "seed")) call get_node_value(root, "seed", seed)
! Check for electron treatment
if (check_for_node(root, "electron_treatment")) then
call get_node_value(root, "electron_treatment", temp_str)
select case (to_lower(temp_str))
case ("led")
electron_treatment = ELECTRON_LED
case ("ttb")
electron_treatment = ELECTRON_TTB
case default
call fatal_error("Unrecognized electron treatment: " // &
trim(temp_str) // ".")
end select
end if
! Check for photon transport
if (check_for_node(root, "photon_transport")) then
call get_node_value(root, "photon_transport", photon_transport)
if (.not. run_CE .and. photon_transport) then
call fatal_error("Photon transport is not currently supported &
&in Multi-group mode")
end if
end if
! Number of bins for logarithmic grid
if (check_for_node(root, "log_grid_bins")) then
call get_node_value(root, "log_grid_bins", n_log_bins)
@ -598,12 +622,18 @@ contains
if (check_for_node(node_cutoff, "weight_avg")) then
call get_node_value(node_cutoff, "weight_avg", weight_survive)
end if
if (check_for_node(node_cutoff, "energy")) then
call get_node_value(node_cutoff, "energy", energy_cutoff(1))
if (check_for_node(node_cutoff, "energy_neutron")) then
call get_node_value(node_cutoff, "energy_neutron", energy_cutoff(1))
end if
if (check_for_node(node_cutoff, "energy_photon")) then
call get_node_value(node_cutoff, "energy_photon", energy_cutoff(2))
end if
if (check_for_node(node_cutoff, "energy_electron")) then
call get_node_value(node_cutoff, "energy_electron", energy_cutoff(3))
end if
if (check_for_node(node_cutoff, "energy_positron")) then
call get_node_value(node_cutoff, "energy_positron", energy_cutoff(4))
end if
end if
! Particle trace
@ -3039,13 +3069,9 @@ contains
! Determine number of bins
select case(temp_str)
case ("energy", "energyout", "mu", "polar", "azimuthal")
if (.not. check_for_node(node_filt, "bins")) then
call fatal_error("Bins not set in filter " // trim(to_str(filter_id)))
end if
n_words = node_word_count(node_filt, "bins")
case ("mesh", "universe", "material", "cell", "distribcell", &
"cellborn", "surface", "delayedgroup")
case ("energy", "energyout", "mu", "polar", "azimuthal", &
"mesh", "universe", "material", "cell", "distribcell", &
"cellborn", "surface", "delayedgroup", "particle")
if (.not. check_for_node(node_filt, "bins")) then
call fatal_error("Bins not set in filter " // trim(to_str(filter_id)))
end if
@ -3099,6 +3125,17 @@ contains
call get_node_array(node_filt, "bins", filt % materials)
end select
case ('particle')
! Allocate and declare the filter type
allocate(ParticleFilter :: f % obj)
select type (filt => f % obj)
type is (ParticleFilter)
! Allocate and store bins
filt % n_bins = n_words
allocate(filt % particles(n_words))
call get_node_array(node_filt, "bins", filt % particles)
end select
case ('universe')
! Allocate and declare the filter type
allocate(UniverseFilter :: f % obj)
@ -3437,6 +3474,8 @@ contains
t % find_filter(FILTER_CELLBORN) = j
type is (MaterialFilter)
t % find_filter(FILTER_MATERIAL) = j
type is (ParticleFilter)
t % find_filter(FILTER_PARTICLE) = j
type is (UniverseFilter)
t % find_filter(FILTER_UNIVERSE) = j
type is (SurfaceFilter)
@ -4071,6 +4110,49 @@ contains
end do
j = j + n_bins
end do
! Check if tally is compatible with particle type
if (photon_transport) then
if (t % find_filter(FILTER_PARTICLE) == 0) then
do j = 1, n_scores
select case (t % score_bins(j))
case (SCORE_INVERSE_VELOCITY)
call fatal_error("Particle filter must be used with photon &
&transport on and inverse velocity score")
case (SCORE_FLUX, SCORE_TOTAL, SCORE_SCATTER, SCORE_NU_SCATTER, &
SCORE_SCATTER_N, SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN, &
SCORE_ABSORPTION, SCORE_FISSION, SCORE_NU_FISSION, &
SCORE_CURRENT, SCORE_FLUX_YN, SCORE_SCATTER_YN, &
SCORE_NU_SCATTER_YN, SCORE_EVENTS, SCORE_DELAYED_NU_FISSION, &
SCORE_PROMPT_NU_FISSION, SCORE_DECAY_RATE)
call warning("Particle filter is not used with photon transport&
& on and " // trim(to_str(t % score_bins(j))) // " score")
end select
end do
else
select type(filt => filters(t % find_filter(FILTER_PARTICLE)) % obj)
type is (ParticleFilter)
do l = 1, filt % n_bins
if (filt % particles(l) == ELECTRON .or. filt % particles(l) == POSITRON) then
t % estimator = ESTIMATOR_ANALOG
end if
end do
end select
end if
else
if (t % find_filter(FILTER_PARTICLE) > 0) then
select type(filt => filters(t % find_filter(FILTER_PARTICLE)) % obj)
type is (ParticleFilter)
do l = 1, filt % n_bins
if (filt % particles(l) /= NEUTRON) then
call warning("Particle filter other than NEUTRON used with &
&photon transport turn off. All tallies for particle &
&type " // trim(to_str(filt % particles(l))) // " will have no scores")
end if
end do
end select
end if
end if
else
call fatal_error("No <scores> specified on tally " &
// trim(to_str(t % id)) // ".")

65
src/nuclide_header.F90
View file

@ -37,12 +37,13 @@ module nuclide_header
end type EnergyGrid
type SumXS
real(8), allocatable :: total(:) ! total cross section
real(8), allocatable :: elastic(:) ! elastic scattering
real(8), allocatable :: fission(:) ! fission
real(8), allocatable :: nu_fission(:) ! neutron production
real(8), allocatable :: absorption(:) ! absorption (MT > 100)
real(8), allocatable :: heating(:) ! heating
real(8), allocatable :: total(:) ! total cross section
real(8), allocatable :: elastic(:) ! elastic scattering
real(8), allocatable :: fission(:) ! fission
real(8), allocatable :: nu_fission(:) ! neutron production
real(8), allocatable :: absorption(:) ! absorption (MT > 100)
real(8), allocatable :: heating(:) ! heating
real(8), allocatable :: photon_prod(:) ! photon production
end type SumXS
type :: Nuclide
@ -91,8 +92,8 @@ module nuclide_header
! array; used at tally-time
! Fission energy release
class(Function1D), allocatable :: fission_q_prompt ! prompt neutrons, gammas
class(Function1D), allocatable :: fission_q_recov ! neutrons, gammas, betas
class(Function1D), allocatable :: fission_q_prompt ! fragments and prompt neutrons, gammas
class(Function1D), allocatable :: fission_q_recov ! fragments, neutrons, gammas, betas
contains
procedure :: clear => nuclide_clear
@ -117,6 +118,7 @@ module nuclide_header
real(8) :: absorption ! microscopic absorption xs
real(8) :: fission ! microscopic fission xs
real(8) :: nu_fission ! microscopic production xs
real(8) :: photon_prod ! microscopic photon production xs
! Information for S(a,b) use
integer :: index_sab ! index in sab_tables (zero means no table)
@ -143,6 +145,7 @@ module nuclide_header
real(8) :: absorption ! macroscopic absorption xs
real(8) :: fission ! macroscopic fission xs
real(8) :: nu_fission ! macroscopic production xs
real(8) :: photon_prod ! macroscopic photon production xs
! Photon cross sections
real(8) :: coherent ! macroscopic coherent xs
@ -442,34 +445,36 @@ contains
if (object_exists(group_id, 'fission_energy_release')) then
fer_group = open_group(group_id, 'fission_energy_release')
! Check to see if this is polynomial or tabulated data
! Q-PROMPT
fer_dset = open_dataset(fer_group, 'q_prompt')
call read_attribute(temp_str, fer_dset, 'type')
if (temp_str == 'Polynomial') then
! Read the prompt Q-value
allocate(Polynomial :: this % fission_q_prompt)
call this % fission_q_prompt % from_hdf5(fer_dset)
call close_dataset(fer_dset)
! Read the recoverable energy Q-value
allocate(Polynomial :: this % fission_q_recov)
fer_dset = open_dataset(fer_group, 'q_recoverable')
call this % fission_q_recov % from_hdf5(fer_dset)
call close_dataset(fer_dset)
else if (temp_str == 'Tabulated1D') then
! Read the prompt Q-value
allocate(Tabulated1D :: this % fission_q_prompt)
call this % fission_q_prompt % from_hdf5(fer_dset)
call close_dataset(fer_dset)
else
call fatal_error('Unrecognized fission prompt energy release format.')
end if
! Read the recoverable energy Q-value
! Q-RECOV
fer_dset = open_dataset(fer_group, 'q_recoverable')
call read_attribute(temp_str, fer_dset, 'type')
if (temp_str == 'Polynomial') then
allocate(Polynomial :: this % fission_q_recov)
call this % fission_q_recov % from_hdf5(fer_dset)
call close_dataset(fer_dset)
else if (temp_str == 'Tabulated1D') then
allocate(Tabulated1D :: this % fission_q_recov)
fer_dset = open_dataset(fer_group, 'q_recoverable')
call this % fission_q_recov % from_hdf5(fer_dset)
call close_dataset(fer_dset)
else
call fatal_error('Unrecognized fission energy release format.')
call fatal_error('Unrecognized fission recoverable energy release format.')
end if
call close_group(fer_group)
end if
@ -481,7 +486,7 @@ contains
subroutine nuclide_create_derived(this)
class(Nuclide), intent(inout) :: this
integer :: i, j, k
integer :: i, j, k, l
integer :: t
integer :: m
integer :: n
@ -501,11 +506,13 @@ contains
allocate(this % sum_xs(i) % fission(n_grid))
allocate(this % sum_xs(i) % nu_fission(n_grid))
allocate(this % sum_xs(i) % absorption(n_grid))
allocate(this % sum_xs(i) % photon_prod(n_grid))
this % sum_xs(i) % total(:) = ZERO
this % sum_xs(i) % elastic(:) = ZERO
this % sum_xs(i) % fission(:) = ZERO
this % sum_xs(i) % nu_fission(:) = ZERO
this % sum_xs(i) % absorption(:) = ZERO
this % sum_xs(i) % photon_prod(:) = ZERO
end do
i_fission = 0
@ -539,6 +546,18 @@ contains
this % sum_xs(t) % total(j:j+n-1) = this % sum_xs(t) % total(j:j+n-1) + &
rx % xs(t) % value
! Calculate nu-photon total cross section
do k = 1, size(rx % products)
if (rx % products(k) % particle == PHOTON) then
do l = 1, n
this % sum_xs(t) % photon_prod(l+j-1) = &
this % sum_xs(t) % photon_prod(l+j-1) + &
rx % xs(t) % value(l) * rx % products(k) % &
yield % evaluate(this % grid(t) % energy(l+j-1))
end do
end if
end do
! Add contribution to absorption cross section
if (is_disappearance(rx % MT)) then
this % sum_xs(t) % absorption(j:j+n-1) = this % sum_xs(t) % &
@ -565,10 +584,6 @@ contains
this % sum_xs(t) % absorption(j:j+n-1) = this % sum_xs(t) % &
absorption(j:j+n-1) + rx % xs(t) % value
! If total fission reaction is present, there's no need to store the
! reaction cross-section since it was copied to this % fission
if (rx % MT == N_FISSION) deallocate(rx % xs(t) % value)
! Keep track of this reaction for easy searching later
if (t == 1) then
i_fission = i_fission + 1

60
src/photon_physics.F90
View file

@ -310,8 +310,17 @@ contains
real(8) :: phi
real(8) :: uvw(3)
! Check for no transitions
if (elm % shells(i_shell) % n_transitions == 0) return
! If no transitions, assume fluorescent photon from captured free electron
if (elm % shells(i_shell) % n_transitions == 0) then
mu = TWO*prn() - ONE
phi = TWO*PI*prn()
uvw(1) = mu
uvw(2) = sqrt(ONE - mu*mu)*cos(phi)
uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
E = elm % shells(i_shell) % binding_energy
call p % create_secondary(uvw, E, PHOTON, run_ce=.true.)
return
end if
! Sample transition
rn = prn()
@ -326,39 +335,36 @@ contains
primary = elm % shells(i_shell) % transition_subshells(1, i_transition)
secondary = elm % shells(i_shell) % transition_subshells(2, i_transition)
if (secondary == 0) then
! Non-radiative trnasition -- Auger/Coster-Kronig effect
! Sample angle isotropically
mu = TWO*prn() - ONE
phi = TWO*PI*prn()
uvw(1) = mu
uvw(2) = sqrt(ONE - mu*mu)*cos(phi)
uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
! TODO: Create electron
! E_electron = transition_energy(i_transition)
! Get the transition energy
E = elm % shells(i_shell) % transition_energy(i_transition)
! Fill secondary (higher) hole first
if (elm % shell_dict % has_key(secondary)) then
i_hole = elm % shell_dict % get_key(secondary)
call atomic_relaxation(p, elm, i_hole)
end if
if (secondary /= 0) then
! Non-radiative transition -- Auger/Coster-Kronig effect
! Create auger electron
call p % create_secondary(uvw, E, ELECTRON, run_ce=.true.)
! Fill hole left by emitted auger electron
i_hole = elm % shell_dict % get_key(secondary)
call atomic_relaxation(p, elm, i_hole)
else
! Radiative transition -- get X-ray energy
E = elm % shells(i_shell) % transition_energy(i_transition)
if (E > ZERO) then
! Sample angle isotropically for X-ray
mu = TWO*prn() - ONE
phi = TWO*PI*prn()
uvw(1) = mu
uvw(2) = sqrt(ONE - mu*mu)*cos(phi)
uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
! Create fluorescent photon
call p % create_secondary(uvw, E, PHOTON, run_ce=.true.)
! Create X-ray
call p % create_secondary(uvw, E, PHOTON, run_ce=.true.)
end if
end if
! Fill primary hole
if (elm % shell_dict % has_key(primary)) then
i_hole = elm % shell_dict % get_key(primary)
call atomic_relaxation(p, elm, i_hole)
end if
! Fill hole created by electron transitioning to the photoelectron hole
i_hole = elm % shell_dict % get_key(primary)
call atomic_relaxation(p, elm, i_hole)
end subroutine atomic_relaxation

247
src/physics.F90
View file

@ -46,15 +46,12 @@ contains
! Sample reaction for the material the particle is in
if (p % type == NEUTRON) then
call sample_neutron_reaction(p)
else
else if (p % type == PHOTON) then
call sample_photon_reaction(p)
end if
! Kill particle if energy falls below cutoff
if (p % E < energy_cutoff(p % type)) then
p % alive = .false.
p % wgt = ZERO
p % last_wgt = ZERO
else if (p % type == ELECTRON) then
call sample_electron_reaction(p)
else if (p % type == POSITRON) then
call sample_positron_reaction(p)
end if
! Display information about collision
@ -113,6 +110,13 @@ contains
end if
end if
! Create secondary photons
if (photon_transport) then
call prn_set_stream(STREAM_PHOTON)
call sample_secondary_photons(p, i_nuclide)
call prn_set_stream(STREAM_TRACKING)
end if
! If survival biasing is being used, the following subroutine adjusts the
! weight of the particle. Otherwise, it checks to see if absorption occurs
@ -159,11 +163,22 @@ contains
real(8) :: cutoff ! sampled total cross section
real(8) :: f ! interpolation factor
real(8) :: xs ! photoionization cross section
real(8) :: r ! random number
real(8) :: prob_after
real(8) :: alpha ! photon energy divided by electron rest mass
real(8) :: alpha_out ! outgoing photon energy over electron rest mass
real(8) :: mu ! scattering cosine
real(8) :: phi ! azimuthal angle
real(8) :: E_electron ! electron energy
real(8) :: uvw(3) ! new direction
real(8) :: rel_vel ! relative velocity of electron
! Kill photon if below energy cutoff
if (p % E < energy_cutoff(PHOTON)) then
p % E = ZERO
p % alive = .false.
return
end if
! Sample element within material
i_element = sample_element(p)
@ -216,12 +231,37 @@ contains
prob = prob + xs
if (prob > cutoff) then
! TODO: Create electron
! E_electron = p % E - elm % shells(i_shell) % binding_energy
E_electron = p % E - elm % shells(i_shell) % binding_energy
! Sample mu using non-relativistic Sauter distribution.
! See Eqns 3.19 and 3.20 in "Implementing a photon physics
! model in Serpent 2" by Toni Kaltiaisenaho
SAMPLE_MU: do
r = prn()
if (FOUR * (ONE - r) * r >= prn()) then
rel_vel = sqrt(E_electron * (E_electron + TWO * MASS_ELECTRON))&
/ (E_electron + MASS_ELECTRON)
mu = (TWO * r + rel_vel - ONE) / &
(TWO * rel_vel * r - rel_vel + ONE)
exit SAMPLE_MU
end if
end do SAMPLE_MU
phi = TWO*PI*prn()
uvw(1) = mu
uvw(2) = sqrt(ONE - mu*mu)*cos(phi)
uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
! Create secondary electron
call p % create_secondary(uvw, E_electron, ELECTRON, run_CE=.true.)
! Allow electrons to fill orbital and produce auger electrons
! and fluorescent photons
call atomic_relaxation(p, elm, i_shell)
p % event_MT = 533 + elm % shells(i_shell) % index_subshell
p % alive = .false.
p % E = ZERO
return
end if
end do
@ -235,21 +275,85 @@ contains
! Sample angle isotropically
mu = TWO*prn() - ONE
phi = TWO*PI*prn()
p % coord(1) % uvw(1) = mu
p % coord(1) % uvw(2) = sqrt(ONE - mu*mu)*cos(phi)
p % coord(1) % uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
uvw(1) = mu
uvw(2) = sqrt(ONE - mu*mu)*cos(phi)
uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
! Set energy
p % E = MASS_ELECTRON
! Compute the kinetic energy of each particle
E_electron = HALF * (p % E - 2 * MASS_ELECTRON)
! Create electron-positron pair traveling in opposite directions
call p % create_secondary( uvw, E_electron, ELECTRON, .true.)
call p % create_secondary(-uvw, E_electron, POSITRON, .true.)
p % event_MT = PAIR_PROD
! Create photon in opposite direction
call p % create_secondary(-p % coord(1) % uvw, MASS_ELECTRON, &
PHOTON, .true.)
p % alive = .false.
p % E = ZERO
end if
end subroutine sample_photon_reaction
!===============================================================================
! SAMPLE_ELECTRON_REACTION terminates the particle and either deposits all
! energy locally (electron_treatment = ELECTRON_LED) or creates secondary
! bremsstrahlung photons from electron deflections with charged particles
! (electron_treatment = ELECTRON_TTB).
!===============================================================================
subroutine sample_electron_reaction(p)
type(Particle), intent(inout) :: p
! TODO: create reaction types
if (electron_treatment == ELECTRON_TTB) then
! TODO: implement thick-target bremsstrahlung model
call fatal_error("Thick-target bremsstrahlung treatment of electrons &
&is not yet implemented.")
end if
p % E = ZERO
p % alive = .false.
end subroutine sample_electron_reaction
!===============================================================================
! SAMPLE_POSITRON_REACTION terminates the particle and either deposits all
! energy locally (electron_treatment = ELECTRON_LED) or creates secondary
! bremsstrahlung photons from electron deflections with charged particles
! (electron_treatment = ELECTRON_TTB). Two annihilation photons of energy
! MASS_ELECTRON (0.511 MeV) are created and travel in opposite directions.
!===============================================================================
subroutine sample_positron_reaction(p)
type(Particle), intent(inout) :: p
real(8) :: mu ! scattering cosine
real(8) :: phi ! azimuthal angle
real(8) :: uvw(3) ! new direction
! TODO: create reaction types
if (electron_treatment == ELECTRON_TTB) then
! TODO: implement thick-target bremsstrahlung model
call fatal_error("Thick-target bremsstrahlung treatment of electrons &
&is not yet implemented.")
end if
! Sample angle isotropically
mu = TWO*prn() - ONE
phi = TWO*PI*prn()
uvw(1) = mu
uvw(2) = sqrt(ONE - mu*mu)*cos(phi)
uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
! Create annihilation photon pair traveling in opposite directions
call p % create_secondary( uvw, MASS_ELECTRON, PHOTON, .true.)
call p % create_secondary(-uvw, MASS_ELECTRON, PHOTON, .true.)
p % E = ZERO
p % alive = .false.
end subroutine sample_positron_reaction
!===============================================================================
! SAMPLE_NUCLIDE
!===============================================================================
@ -421,6 +525,66 @@ contains
end subroutine sample_fission
!===============================================================================
! SAMPLE_PHOTON_PRODUCT
!===============================================================================
subroutine sample_photon_product(i_nuclide, E, i_reaction, i_product)
integer, intent(in) :: i_nuclide ! index in nuclides array
real(8), intent(in) :: E ! energy of neutron
integer, intent(out) :: i_reaction ! index in nuc % reactions array
integer, intent(out) :: i_product ! index in nuc % reactions array
integer :: i_grid
integer :: i_temp
integer :: threshold
integer :: last_valid_reaction
integer :: last_valid_product
real(8) :: f
real(8) :: prob
real(8) :: cutoff
real(8) :: yield
! Get pointer to nuclide
associate (nuc => nuclides(i_nuclide))
! Get grid index and interpolation factor and sample proton production cdf
i_temp = micro_xs(i_nuclide) % index_temp
i_grid = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
cutoff = prn() * micro_xs(i_nuclide) % photon_prod
prob = ZERO
! Loop through each reaction type
REACTION_LOOP: do i_reaction = 1, size(nuc % reactions)
associate (rx => nuc % reactions(i_reaction))
do i_product = 1, size(rx % products)
if (rx % products(i_product) % particle == PHOTON) then
threshold = rx % xs(i_temp) % threshold
! if energy is below threshold for this reaction, skip it
if (i_grid < threshold) cycle
! add to cumulative probability
yield = rx % products(i_product) % yield % evaluate(E)
prob = prob + ((ONE - f) * rx % xs(i_temp) % value(i_grid - threshold + 1) &
+ f*(rx % xs(i_temp) % value(i_grid - threshold + 2))) * yield
if (prob > cutoff) return
last_valid_reaction = i_reaction
last_valid_product = i_product
end if
end do
end associate
end do REACTION_LOOP
end associate
i_reaction = last_valid_reaction
i_product = last_valid_product
end subroutine sample_photon_product
!===============================================================================
! ABSORPTION
!===============================================================================
@ -1498,4 +1662,49 @@ contains
end subroutine inelastic_scatter
!===============================================================================
! SAMPLE_SECONDARY_PHOTONS
!===============================================================================
subroutine sample_secondary_photons(p, i_nuclide)
type(Particle), intent(inout) :: p
integer, intent(in) :: i_nuclide
integer :: i_reaction ! index in nuc % reactions array
integer :: i_product ! index in nuc % reactions % products array
real(8) :: nu_t
real(8) :: mu
real(8) :: E
real(8) :: uvw(3)
integer :: nu
integer :: i
! Sample the number of photons produced
nu_t = micro_xs(i_nuclide) % photon_prod / micro_xs(i_nuclide) % total
if (prn() > nu_t - int(nu_t)) then
nu = int(nu_t)
else
nu = int(nu_t) + 1
end if
! Sample each secondary photon
do i = 1, nu
! Sample the reaction and product
call sample_photon_product(i_nuclide, p % E, i_reaction, i_product)
! Sample the outgoing energy and angle
call nuclides(i_nuclide) % reactions(i_reaction) % products(i_product) &
% sample(p % E, E, mu)
! Sample the new direction
uvw = rotate_angle(p % coord(1) % uvw, mu)
! Create the secondary photon
call p % create_secondary(uvw, E, PHOTON, run_CE=.true.)
end do
end subroutine sample_secondary_photons
end module physics

151
src/tally.F90
View file

@ -99,6 +99,7 @@ contains
real(8) :: f ! interpolation factor
real(8) :: score ! analog tally score
real(8) :: E ! particle energy
real(8) :: xs ! cross section
i = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
@ -117,6 +118,7 @@ contains
case (SCORE_FLUX, SCORE_FLUX_YN)
if (t % estimator == ESTIMATOR_ANALOG) then
! All events score to a flux bin. We actually use a collision
! estimator in place of an analog one since there is no way to count
@ -128,7 +130,12 @@ contains
else
score = p % last_wgt
end if
score = score / material_xs % total * flux
if (p % type == NEUTRON .or. p % type == PHOTON) then
score = score / material_xs % total * flux
else
score = ZERO
end if
else
! For flux, we need no cross section
@ -159,7 +166,6 @@ contains
case (SCORE_INVERSE_VELOCITY)
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
E = p % E
@ -339,6 +345,9 @@ contains
case (SCORE_FISSION)
if (material_xs % absorption == ZERO) cycle SCORE_LOOP
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
@ -370,11 +379,15 @@ contains
case (SCORE_NU_FISSION)
if (material_xs % absorption == ZERO) cycle SCORE_LOOP
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing .or. p % fission) then
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
! Normally, we only need to make contributions to one scoring
! bin. However, in the case of fission, since multiple fission
! neutrons were emitted with different energies, multiple
! outgoing energy bins may have been scored to. The following
! logic treats this special case and results to multiple bins
@ -413,6 +426,9 @@ contains
case (SCORE_PROMPT_NU_FISSION)
if (material_xs % absorption == ZERO) cycle SCORE_LOOP
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
E = p % E
@ -484,6 +500,8 @@ contains
case (SCORE_DELAYED_NU_FISSION)
if (material_xs % absorption == ZERO) cycle SCORE_LOOP
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
E = p % E
@ -683,6 +701,8 @@ contains
case (SCORE_DECAY_RATE)
if (material_xs % absorption == ZERO) cycle SCORE_LOOP
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
E = p % E
@ -969,6 +989,9 @@ contains
end if
case (SCORE_KAPPA_FISSION)
if (material_xs % absorption == ZERO) cycle SCORE_LOOP
! 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))
@ -1053,7 +1076,10 @@ contains
end if
end if
case (SCORE_FISS_Q_PROMPT)
case (SCORE_FISS_Q_PROMPT, SCORE_FISS_Q_RECOV)
if (material_xs % absorption == ZERO) cycle SCORE_LOOP
score = ZERO
if (t % estimator == ESTIMATOR_ANALOG) then
@ -1064,10 +1090,15 @@ contains
associate (nuc => nuclides(p % event_nuclide))
if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
allocated(nuc % fission_q_prompt)) then
score = p % absorb_wgt &
* nuc % fission_q_prompt % evaluate(p % last_E) &
if (score_bin == SCORE_FISS_Q_PROMPT) then
xs = nuc % fission_q_prompt % evaluate(p % last_E)
else if (score_bin == SCORE_FISS_Q_RECOV) then
xs = nuc % fission_q_recov % evaluate(p % last_E)
end if
score = p % absorb_wgt * xs * flux &
* micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption * flux
/ micro_xs(p % event_nuclide) % absorption
end if
end associate
else
@ -1076,72 +1107,19 @@ contains
! 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
associate (nuc => nuclides(p % event_nuclide))
if (allocated(nuc % fission_q_prompt)) then
score = p % last_wgt &
* nuc % fission_q_prompt % evaluate(p % last_E) &
* micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption * flux
end if
end associate
end if
else
if (t % estimator == ESTIMATOR_COLLISION) then
E = p % last_E
else
E = p % E
end if
if (i_nuclide > 0) then
if (allocated(nuclides(i_nuclide) % fission_q_prompt)) then
score = micro_xs(i_nuclide) % fission * atom_density * flux &
* nuclides(i_nuclide) % fission_q_prompt % evaluate(E)
end if
else
if (p % material /= MATERIAL_VOID) then
do l = 1, materials(p % material) % n_nuclides
atom_density_ = materials(p % material) % atom_density(l)
i_nuc = materials(p % material) % nuclide(l)
if (allocated(nuclides(i_nuc) % fission_q_prompt)) then
score = score + micro_xs(i_nuc) % fission * atom_density_ &
* flux &
* nuclides(i_nuc) % fission_q_prompt % evaluate(E)
end if
end do
end if
end if
end if
case (SCORE_FISS_Q_RECOV)
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 scaled by Q-value
associate (nuc => nuclides(p % event_nuclide))
if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
allocated(nuc % fission_q_recov)) then
score = p % absorb_wgt &
* nuc % fission_q_recov % evaluate(p % last_E) &
allocated(nuc % fission_q_prompt)) then
if (score_bin == SCORE_FISS_Q_PROMPT) then
xs = nuc % fission_q_prompt % evaluate(p % last_E)
else if (score_bin == SCORE_FISS_Q_RECOV) then
xs = nuc % fission_q_recov % evaluate(p % last_E)
end if
score = p % last_wgt * xs * flux &
* micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption * flux
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
associate (nuc => nuclides(p % event_nuclide))
if (allocated(nuc % fission_q_recov)) then
score = p % last_wgt &
* nuc % fission_q_recov % evaluate(p % last_E) &
* micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption * flux
/ micro_xs(p % event_nuclide) % absorption
end if
end associate
end if
@ -1154,20 +1132,37 @@ contains
end if
if (i_nuclide > 0) then
if (allocated(nuclides(i_nuclide) % fission_q_recov)) then
score = micro_xs(i_nuclide) % fission * atom_density * flux &
* nuclides(i_nuclide) % fission_q_recov % evaluate(E)
end if
associate (nuc => nuclides(i_nuclide))
if (allocated(nuc % fission_q_prompt)) then
if (score_bin == SCORE_FISS_Q_PROMPT) then
xs = nuc % fission_q_prompt % evaluate(E)
else if (score_bin == SCORE_FISS_Q_RECOV) then
xs = nuc % fission_q_recov % evaluate(E)
end if
score = micro_xs(i_nuclide) % fission * atom_density * flux * xs
end if
end associate
else
if (p % material /= MATERIAL_VOID) then
do l = 1, materials(p % material) % n_nuclides
atom_density_ = materials(p % material) % atom_density(l)
i_nuc = materials(p % material) % nuclide(l)
if (allocated(nuclides(i_nuc) % fission_q_recov)) then
score = score + micro_xs(i_nuc) % fission * atom_density_ &
* flux &
* nuclides(i_nuc) % fission_q_recov % evaluate(E)
end if
associate (nuc => nuclides(i_nuc))
if (allocated(nuc % fission_q_prompt)) then
if (score_bin == SCORE_FISS_Q_PROMPT) then
xs = nuc % fission_q_prompt % evaluate(E)
else if (score_bin == SCORE_FISS_Q_RECOV) then
xs = nuc % fission_q_recov % evaluate(E)
end if
score = score + micro_xs(i_nuc) % fission * atom_density_ &
* flux * xs
end if
end associate
end do
end if
end if

71
src/tally_filter.F90
View file

@ -58,6 +58,19 @@ module tally_filter
procedure :: initialize => initialize_material
end type MaterialFilter
!===============================================================================
! PARTICLE specifies which particle tally events reside in.
!===============================================================================
type, extends(TallyFilter) :: ParticleFilter
integer, allocatable :: particles(:)
type(DictIntInt) :: map
contains
procedure :: get_next_bin => get_next_bin_particle
procedure :: to_statepoint => to_statepoint_particle
procedure :: text_label => text_label_particle
procedure :: initialize => initialize_particle
end type ParticleFilter
!===============================================================================
! CELLFILTER specifies which geometric cells tally events reside in.
!===============================================================================
@ -313,6 +326,12 @@ contains
! Compute the length of the entire track.
total_distance = sqrt(sum((xyz1 - xyz0)**2))
! Check if particle has moved
if (total_distance == ZERO) then
next_bin = current_bin
return
end if
if (current_bin == NO_BIN_FOUND) then
! We are looking for the first valid mesh bin. Check to see if the
! particle starts inside the mesh.
@ -632,6 +651,58 @@ contains
label = "Material " // to_str(materials(this % materials(bin)) % id)
end function text_label_material
!===============================================================================
! ParticleFilter methods
!===============================================================================
subroutine get_next_bin_particle(this, p, estimator, current_bin, next_bin, &
weight)
class(ParticleFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
integer, value, intent(in) :: current_bin
integer, intent(out) :: next_bin
real(8), intent(out) :: weight
integer :: i
weight = ERROR_REAL
next_bin = NO_BIN_FOUND
if (current_bin == NO_BIN_FOUND) then
do i = 1, this % n_bins
if (this % particles(i) == p % type) then
next_bin = i
weight = ONE
end if
end do
end if
end subroutine get_next_bin_particle
subroutine to_statepoint_particle(this, filter_group)
class(ParticleFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
integer :: i
call write_dataset(filter_group, "type", "particle")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % particles)
end subroutine to_statepoint_particle
subroutine initialize_particle(this)
class(ParticleFilter), intent(inout) :: this
end subroutine initialize_particle
function text_label_particle(this, bin) result(label)
class(ParticleFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
label = "Particle " // to_str(this % particles(bin))
end function text_label_particle
!===============================================================================
! CellFilter methods
!===============================================================================

19
src/tracking.F90
View file

@ -11,7 +11,7 @@ module tracking
use particle_header, only: LocalCoord, Particle
use physics, only: collision
use physics_mg, only: collision_mg
use random_lcg, only: prn
use random_lcg, only: prn, prn_set_stream
use string, only: to_str
use tally, only: score_analog_tally, score_tracklength_tally, &
score_collision_tally, score_surface_current, &
@ -69,6 +69,14 @@ contains
if (active_tallies % size() > 0) call zero_flux_derivs()
EVENT_LOOP: do
! Set the random number stream
if (p % type == NEUTRON) then
call prn_set_stream(STREAM_TRACKING)
else
call prn_set_stream(STREAM_PHOTON)
end if
! If the cell hasn't been determined based on the particle's location,
! initiate a search for the current cell. This generally happens at the
! beginning of the history and again for any secondary particles
@ -114,7 +122,9 @@ contains
lattice_translation, next_level)
! Sample a distance to collision
if (material_xs % total == ZERO) then
if (p % type == ELECTRON .or. p % type == POSITRON) then
d_collision = ZERO
else if (material_xs % total == ZERO) then
d_collision = INFINITY
else
d_collision = -log(prn()) / material_xs % total
@ -133,9 +143,8 @@ contains
call score_tracklength_tally(p, distance)
end if
! Score track-length estimate of k-eff
if (run_mode == MODE_EIGENVALUE) then
if (run_mode == MODE_EIGENVALUE .and. p % type == NEUTRON) then
global_tally_tracklength = global_tally_tracklength + p % wgt * &
distance * material_xs % nu_fission
end if
@ -166,7 +175,7 @@ contains
! PARTICLE HAS COLLISION
! Score collision estimate of keff
if (run_mode == MODE_EIGENVALUE) then
if (run_mode == MODE_EIGENVALUE .and. p % type == NEUTRON) then
global_tally_collision = global_tally_collision + p % wgt * &
material_xs % nu_fission / material_xs % total
end if