diff --git a/openmc/data/library.py b/openmc/data/library.py index 2f33743a6e..baf8f6ae0b 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -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 diff --git a/openmc/filter.py b/openmc/filter.py index 90f749668d..d971163ee0 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -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. diff --git a/openmc/settings.py b/openmc/settings.py index 658dd8642c..84f77312d9 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -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) diff --git a/scripts/openmc-get-nndc-data b/scripts/openmc-get-nndc-data index f1da2241ab..8afd98cf70 100755 --- a/scripts/openmc-get-nndc-data +++ b/scripts/openmc-get-nndc-data @@ -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']) + diff --git a/scripts/openmc-get-photo-endf71 b/scripts/openmc-get-photo-endf71 new file mode 100755 index 0000000000..9920b27e28 --- /dev/null +++ b/scripts/openmc-get-photo-endf71 @@ -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) diff --git a/src/constants.F90 b/src/constants.F90 index 53253e23db..16d266a5a1 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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 diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 918b06ad83..37a3768455 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -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 diff --git a/src/energy_distribution.F90 b/src/energy_distribution.F90 index 770da617cc..abfdf6f2f5 100644 --- a/src/energy_distribution.F90 +++ b/src/energy_distribution.F90 @@ -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 diff --git a/src/global.F90 b/src/global.F90 index d9da251bad..30d824a6d2 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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 ! ============================================================================ diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 0dc1a90a83..5273a21807 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 specified on tally " & // trim(to_str(t % id)) // ".") diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 9cc2bc7e20..ab937809ee 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -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 diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 7b7244fe70..7d4809803d 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -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 diff --git a/src/physics.F90 b/src/physics.F90 index 7582190834..6923a5386a 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 diff --git a/src/tally.F90 b/src/tally.F90 index 14940551b7..3d7ee2cb80 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -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 diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index dfb89d4330..a010a61d5d 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -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 !=============================================================================== diff --git a/src/tracking.F90 b/src/tracking.F90 index 7650f598ab..c4baa2c012 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -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