From 8da0636ab827b1d409425080aa81ef215b738be6 Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Thu, 13 Jul 2017 09:32:02 -0400 Subject: [PATCH] added particle filter and addressed PR comments --- openmc/filter.py | 40 ++- openmc/settings.py | 26 +- scripts/openmc-get-photo-endf71 | 58 ++-- src/constants.F90 | 32 +- src/cross_section.F90 | 13 +- src/global.F90 | 3 +- src/input_xml.F90 | 118 ++++++-- src/nuclide_header.F90 | 232 +++++++-------- src/output.F90 | 20 +- src/physics.F90 | 66 +++-- src/tally.F90 | 506 +++++++++++++++----------------- src/tally_filter.F90 | 65 ++++ src/tracking.F90 | 10 +- 13 files changed, 687 insertions(+), 502 deletions(-) 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-photo-endf71 b/scripts/openmc-get-photo-endf71 index b2ee3bba2a..7fd06a0a24 100755 --- a/scripts/openmc-get-photo-endf71 +++ b/scripts/openmc-get-photo-endf71 @@ -2,6 +2,7 @@ from __future__ import print_function import os +import sys import shutil import zipfile import requests @@ -31,57 +32,48 @@ parser.add_argument('-c', '--cross-sections-file', help='cross_sections.xml file to append libraries to') args = parser.parse_args() -base_url = 'http://www-nds.iaea.org/public/download-endf/ENDF-B-VII.1/' +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 IAEA SITE AND GENERATE HDF5 LIBRARY - -# Make photo and ard directories -if not os.path.exists('photo'): - os.mkdir('photo') +# DOWNLOAD FILES FROM NNDC SITE if not os.path.exists('photo_hdf5'): os.mkdir('photo_hdf5') -if not os.path.exists('ard'): - os.mkdir('ard') - library = openmc.data.DataLibrary() -for z in range(1,101): +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)) - # Download photo files - if z < 100: - filename = 'photo/photo_{:02}00_{}-{}-0'.format(z, z, element) - else: - filename = 'photo/photo_{}20_{}-{}-0'.format(z-1, z, element) + # Load files + filename = 'photoat/photoat-{:03}_{}_000.endf'.format(z, element) + photo_file = 'photoat/' + element + '.endf' + shutil.move(filename, photo_file) - url = base_url + filename + '.zip' - r = requests.get(url, stream=True) - zipfile.ZipFile(BytesIO(r.content)).extractall(path='photo') - photo_file = 'photo/' + element + '.dat' - shutil.move(filename + '.dat', photo_file) - - # Download ard files - if z < 100: - filename = 'ard/ard_{:02}00_{}-{}-0'.format(z, z, element) - else: - filename = 'ard/ard_{}20_{}-{}-0'.format(z-1, z, element) - - url = base_url + filename + '.zip' - r = requests.get(url, stream=True) - zipfile.ZipFile(BytesIO(r.content)).extractall(path='ard') - ard_file = 'ard/' + element + '.dat' - shutil.move(filename + '.dat', ard_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, ard_file) + f = openmc.data.IncidentPhoton.from_endf(photo_file, atom_file) f.export_to_hdf5(hdf5_file) library.register_file(hdf5_file) diff --git a/src/constants.F90 b/src/constants.F90 index 907a3dbcc8..c04ecf1b2b 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -318,7 +318,7 @@ module constants EVENT_ABSORB = 2 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 30 + integer, parameter :: N_SCORE_TYPES = 36 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -341,15 +341,21 @@ 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_HEATING = -22, & ! prompt fission Q-value SCORE_FISS_Q_RECOV = -23, & ! recoverable fission Q-value - SCORE_FISS_Q_FRAGMENTS = -24, & ! fragment fission Q-value - SCORE_FISS_Q_BETAS = -25, & ! beta fission Q-value - SCORE_Q_ELASTIC = -26, & ! elastic scatter Q-value - SCORE_Q_PHOTONS = -27, & ! photon Q-value below threshold - SCORE_Q_ELECTRONS = -28, & ! electron Q-value - SCORE_Q_POSITRONS = -29, & ! positron Q-value - SCORE_DECAY_RATE = -30 ! delayed neutron precursor decay rate + SCORE_FISS_Q_PROMPT = -24, & ! recoverable fission Q-value + SCORE_FISS_Q_PROMPT_NEUTRONS = -25, & ! recoverable fission Q-value + SCORE_FISS_Q_DELAYED_NEUTRONS = -26, & ! recoverable fission Q-value + SCORE_FISS_Q_FRAGMENTS = -27, & ! recoverable fission Q-value + SCORE_FISS_Q_BETAS = -28, & ! recoverable fission Q-value + SCORE_FISS_Q_PROMPT_PHOTONS = -29, & ! recoverable fission Q-value + SCORE_FISS_Q_DELAYED_PHOTONS = -30, & ! recoverable fission Q-value + SCORE_FISS_Q_NEUTRINOS = -31, & ! recoverable fission Q-value + SCORE_Q_PHOTONS = -32, & ! recoverable fission Q-value + SCORE_Q_ELECTRONS = -33, & ! recoverable fission Q-value + SCORE_Q_POSITRONS = -34, & ! recoverable fission Q-value + SCORE_Q_ELASTIC = -35, & ! recoverable fission Q-value + SCORE_DECAY_RATE = -36 ! delayed neutron precursor decay rate ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 @@ -372,7 +378,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, & @@ -387,7 +393,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 :: & @@ -431,12 +438,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 diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 817a874445..325f1ed841 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -263,7 +263,7 @@ contains ! Initialize nuclide cross-sections to zero micro_xs(i_nuclide) % fission = ZERO micro_xs(i_nuclide) % nu_fission = ZERO - micro_xs(i_nuclide) % nu_photon_total = 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) & @@ -277,9 +277,9 @@ contains micro_xs(i_nuclide) % absorption = (ONE - f) * xs % absorption( & i_grid) + f * xs % absorption(i_grid + 1) - ! Calculate microscopic nuclide nu-photon total cross section - micro_xs(i_nuclide) % nu_photon_total = (ONE - f) * xs % & - nu_photon_total(i_grid) + f * xs % nu_photon_total(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 @@ -573,11 +573,6 @@ contains micro_xs(i_nuclide) % fission = fission micro_xs(i_nuclide) % total = elastic + inelastic + capture + fission - ! Set the nu-photon production cross section - i_grid = int(log(E/energy_min_neutron)/log_spacing) - micro_xs(i_nuclide) % nu_photon_total = & - nuc % compute_nu_photon_total(E, i_temp, i_grid) - ! Determine nu-fission cross section if (nuc % fissionable) then micro_xs(i_nuclide) % nu_fission = nuc % nu(E, EMISSION_TOTAL) * & diff --git a/src/global.F90 b/src/global.F90 index 680fed2cb9..28ff57cc60 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -59,6 +59,7 @@ module global type(DictIntInt) :: lattice_dict type(DictIntInt) :: surface_dict type(DictIntInt) :: material_dict + type(DictIntInt) :: particle_dict type(DictIntInt) :: mesh_dict type(DictIntInt) :: filter_dict type(DictIntInt) :: tally_dict @@ -120,7 +121,7 @@ module global integer :: n_log_bins ! number of bins for logarithmic grid real(8) :: log_spacing ! spacing on logarithmic grid - logical :: photon_transport = .true. + logical :: photon_transport = .false. integer :: electron_treatment = ELECTRON_LED ! ============================================================================ diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 59359ef115..34b14cddec 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -271,9 +271,9 @@ contains 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") + case ("led") electron_treatment = ELECTRON_LED - case ("TTB") + case ("ttb") electron_treatment = ELECTRON_TTB case default call fatal_error("Unrecognized electron treatment: " // & @@ -283,16 +283,12 @@ contains ! Check for photon transport if (check_for_node(root, "photon_transport")) then - call get_node_value(root, "photon_transport", temp_str) - select case (to_lower(temp_str)) - case ("true") - photon_transport = .true. - case ("false") - photon_transport = .false. - case default - call fatal_error("Unrecognized photon transport: " // & - trim(temp_str) // ".") - end select + call get_node_value(root, "photon_transport", photon_transport) + + if (.not. run_CE) then + call fatal_error("Photon transport is not currently supported & + &in Multi-group mode") + end if end if ! Number of bins for logarithmic grid @@ -626,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 @@ -3067,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 @@ -3127,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) @@ -3465,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) @@ -3872,22 +3883,32 @@ contains t % score_bins(j) = SCORE_FISS_Q_PROMPT case ('fission-q-recoverable') t % score_bins(j) = SCORE_FISS_Q_RECOV + case ('fission-q-prompt-neutrons') + t % score_bins(j) = SCORE_FISS_Q_PROMPT_NEUTRONS + case ('fission-q-delayed-neutrons') + t % score_bins(j) = SCORE_FISS_Q_DELAYED_NEUTRONS case ('fission-q-fragments') t % score_bins(j) = SCORE_FISS_Q_FRAGMENTS case ('fission-q-betas') t % score_bins(j) = SCORE_FISS_Q_BETAS - case ('q-elastic') - t % score_bins(j) = SCORE_Q_ELASTIC - t % estimator = ESTIMATOR_ANALOG - case ('q-photons') - t % score_bins(j) = SCORE_Q_PHOTONS - t % estimator = ESTIMATOR_ANALOG + case ('fission-q-prompt-photons') + t % score_bins(j) = SCORE_FISS_Q_PROMPT_PHOTONS + case ('fission-q-delayed-photons') + t % score_bins(j) = SCORE_FISS_Q_DELAYED_PHOTONS + case ('fission-q-neutrinos') + t % score_bins(j) = SCORE_FISS_Q_NEUTRINOS case ('q-electrons') t % score_bins(j) = SCORE_Q_ELECTRONS t % estimator = ESTIMATOR_ANALOG case ('q-positrons') t % score_bins(j) = SCORE_Q_POSITRONS t % estimator = ESTIMATOR_ANALOG + case ('q-elastic') + t % score_bins(j) = SCORE_Q_ELASTIC + t % estimator = ESTIMATOR_ANALOG + case ('heating') + t % score_bins(j) = SCORE_HEATING + t % estimator = ESTIMATOR_ANALOG case ('current') t % score_bins(j) = SCORE_CURRENT t % type = TALLY_SURFACE_CURRENT @@ -4115,6 +4136,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 % bins % data(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 % bins % data(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 716c482067..84de722a7d 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -5,7 +5,7 @@ module nuclide_header use hdf5, only: HID_T, HSIZE_T, SIZE_T - use algorithm, only: sort, find, binary_search + use algorithm, only: sort, find use constants use dict_header, only: DictIntInt use endf, only: reaction_name, is_fission, is_disappearance @@ -43,7 +43,7 @@ module nuclide_header real(8), allocatable :: nu_fission(:) ! neutron production real(8), allocatable :: absorption(:) ! absorption (MT > 100) real(8), allocatable :: heating(:) ! heating - real(8), allocatable :: nu_photon_total(:) ! photon production + real(8), allocatable :: photon_prod(:) ! photon production end type SumXS type :: Nuclide @@ -96,13 +96,17 @@ module nuclide_header class(Function1D), allocatable :: fission_q_recov ! fragments, neutrons, gammas, betas class(Function1D), allocatable :: fission_q_fragments ! fragments class(Function1D), allocatable :: fission_q_betas ! betas + class(Function1D), allocatable :: fission_q_neutrinos ! betas + class(Function1D), allocatable :: fission_q_delayed_neutrons ! betas + class(Function1D), allocatable :: fission_q_prompt_neutrons ! betas + class(Function1D), allocatable :: fission_q_delayed_photons ! betas + class(Function1D), allocatable :: fission_q_prompt_photons ! betas contains procedure :: clear => nuclide_clear procedure :: from_hdf5 => nuclide_from_hdf5 procedure :: init_grid => nuclide_init_grid procedure :: nu => nuclide_nu - procedure :: compute_nu_photon_total => compute_nuclide_nu_photon_total procedure, private :: create_derived => nuclide_create_derived end type Nuclide @@ -121,7 +125,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) :: nu_photon_total ! microscopic photon 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) @@ -148,7 +152,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) :: nu_photon_total ! macroscopic photon production xs + real(8) :: photon_prod ! macroscopic photon production xs ! Photon cross sections real(8) :: coherent ! macroscopic coherent xs @@ -448,79 +452,141 @@ 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 - allocate(Tabulated1D :: this % fission_q_recov) - fer_dset = open_dataset(fer_group, 'q_recoverable') + ! 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) + call this % fission_q_recov % from_hdf5(fer_dset) + call close_dataset(fer_dset) + else + call fatal_error('Unrecognized fission recoverable energy release format.') + end if + + ! Q-FRAGMENTS + fer_dset = open_dataset(fer_group, 'fragments') + call read_attribute(temp_str, fer_dset, 'type') + if (temp_str == 'Polynomial') then + allocate(Polynomial :: this % fission_q_fragments) + call this % fission_q_fragments % from_hdf5(fer_dset) + call close_dataset(fer_dset) + else if (temp_str == 'Tabulated1D') then + allocate(Tabulated1D :: this % fission_q_fragments) + call this % fission_q_fragments % from_hdf5(fer_dset) + call close_dataset(fer_dset) + else + call fatal_error('Unrecognized fission fragments energy release format.') + end if + + ! Q-BETAS + fer_dset = open_dataset(fer_group, 'betas') + call read_attribute(temp_str, fer_dset, 'type') + if (temp_str == 'Polynomial') then + allocate(Polynomial :: this % fission_q_betas) + call this % fission_q_betas % from_hdf5(fer_dset) + call close_dataset(fer_dset) + else if (temp_str == 'Tabulated1D') then + allocate(Tabulated1D :: this % fission_q_betas) + call this % fission_q_betas % from_hdf5(fer_dset) + call close_dataset(fer_dset) + else + call fatal_error('Unrecognized fission betas energy release format.') + end if + + ! Q-NEUTRINOS + fer_dset = open_dataset(fer_group, 'neutrinos') + call read_attribute(temp_str, fer_dset, 'type') + if (temp_str == 'Polynomial') then + allocate(Polynomial :: this % fission_q_neutrinos) + call this % fission_q_neutrinos % from_hdf5(fer_dset) + call close_dataset(fer_dset) + else if (temp_str == 'Tabulated1D') then + allocate(Tabulated1D :: this % fission_q_neutrinos) + call this % fission_q_neutrinos % from_hdf5(fer_dset) + call close_dataset(fer_dset) + else + call fatal_error('Unrecognized fission neutrinos energy release format.') + end if + + ! Q-DELAYED-NEUTRONS + fer_dset = open_dataset(fer_group, 'delayed_neutrons') + call read_attribute(temp_str, fer_dset, 'type') + if (temp_str == 'Polynomial') then + allocate(Polynomial :: this % fission_q_delayed_neutrons) + call this % fission_q_delayed_neutrons % from_hdf5(fer_dset) + call close_dataset(fer_dset) + else if (temp_str == 'Tabulated1D') then + allocate(Tabulated1D :: this % fission_q_delayed_neutrons) + call this % fission_q_delayed_neutrons % from_hdf5(fer_dset) + call close_dataset(fer_dset) + else + call fatal_error('Unrecognized fission delayed neutron energy release format.') + end if + + ! Q-PROMPT-NEUTRONS + fer_dset = open_dataset(fer_group, 'prompt_neutrons') + call read_attribute(temp_str, fer_dset, 'type') + if (temp_str == 'Polynomial') then + allocate(Polynomial :: this % fission_q_prompt_neutrons) + call this % fission_q_prompt_neutrons % from_hdf5(fer_dset) + call close_dataset(fer_dset) + else if (temp_str == 'Tabulated1D') then + allocate(Tabulated1D :: this % fission_q_prompt_neutrons) + call this % fission_q_prompt_neutrons % from_hdf5(fer_dset) + call close_dataset(fer_dset) else call fatal_error('Unrecognized fission energy release format.') end if - fer_dset = open_dataset(fer_group, 'fragments') + ! Q-DELAYED-PHOTONS + fer_dset = open_dataset(fer_group, 'delayed_photons') call read_attribute(temp_str, fer_dset, 'type') if (temp_str == 'Polynomial') then - - ! Read the fragment energy Q-value - allocate(Polynomial :: this % fission_q_fragments) - fer_dset = open_dataset(fer_group, 'fragments') - call this % fission_q_fragments % from_hdf5(fer_dset) + allocate(Polynomial :: this % fission_q_delayed_photons) + call this % fission_q_delayed_photons % from_hdf5(fer_dset) call close_dataset(fer_dset) - else if (temp_str == 'Tabulated1D') then - - ! Read the fragment energy Q-value - allocate(Tabulated1D :: this % fission_q_fragments) - fer_dset = open_dataset(fer_group, 'fragments') - call this % fission_q_fragments % from_hdf5(fer_dset) + allocate(Tabulated1D :: this % fission_q_delayed_photons) + call this % fission_q_delayed_photons % from_hdf5(fer_dset) call close_dataset(fer_dset) else - call fatal_error('Unrecognized fission fragment energy release format.') + call fatal_error('Unrecognized fission delayed photon energy release format.') end if - fer_dset = open_dataset(fer_group, 'betas') + ! Q-PROMPT-PHOTONS + fer_dset = open_dataset(fer_group, 'prompt_photons') call read_attribute(temp_str, fer_dset, 'type') if (temp_str == 'Polynomial') then - - ! Read the beta energy Q-value - allocate(Polynomial :: this % fission_q_betas) - fer_dset = open_dataset(fer_group, 'betas') - call this % fission_q_betas % from_hdf5(fer_dset) + allocate(Polynomial :: this % fission_q_prompt_photons) + call this % fission_q_prompt_photons % from_hdf5(fer_dset) call close_dataset(fer_dset) - else if (temp_str == 'Tabulated1D') then - - ! Read the beta energy Q-value - allocate(Tabulated1D :: this % fission_q_betas) - fer_dset = open_dataset(fer_group, 'betas') - call this % fission_q_betas % from_hdf5(fer_dset) + allocate(Tabulated1D :: this % fission_q_prompt_photons) + call this % fission_q_prompt_photons % from_hdf5(fer_dset) call close_dataset(fer_dset) else - call fatal_error('Unrecognized beta energy release format.') + call fatal_error('Unrecognized fission prompt photon energy release format.') end if + call close_group(fer_group) end if @@ -552,13 +618,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) % nu_photon_total(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) % nu_photon_total(:) = ZERO + this % sum_xs(i) % photon_prod(:) = ZERO end do i_fission = 0 @@ -596,8 +662,8 @@ contains do k = 1, size(rx % products) if (rx % products(k) % particle == PHOTON) then do l = 1, n - this % sum_xs(t) % nu_photon_total(l+j-1) = & - this % sum_xs(t) % nu_photon_total(l+j-1) + & + 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 @@ -751,72 +817,6 @@ contains end function nuclide_nu -!=============================================================================== -! COMPUTE_NUCLIDE_NU_PHOTON_TOTAL is an interface to compute the number of -! photons produced -!=============================================================================== - - pure function compute_nuclide_nu_photon_total(this, E, t, i_log_union) result(nu_photon_total) - class(Nuclide), intent(in) :: this - real(8), intent(in) :: E - integer, intent(in) :: t - integer, intent(in) :: i_log_union - real(8) :: rx_xs - real(8) :: nu_photon_total - real(8) :: f - integer :: m, j, k - integer :: i_grid, i_low, i_high - - associate (grid => this % grid(t), xs => this % sum_xs(t)) - ! Determine the energy grid index using a logarithmic mapping to - ! reduce the energy range over which a binary search needs to be - ! performed - - if (E < grid % energy(1)) then - i_grid = 1 - elseif (E > grid % energy(size(grid % energy))) then - i_grid = size(grid % energy) - 1 - else - ! Determine bounding indices based on which equal log-spaced - ! interval the energy is in - i_low = grid % grid_index(i_log_union) - i_high = grid % grid_index(i_log_union + 1) + 1 - - ! Perform binary search over reduced range - i_grid = binary_search(grid % energy(i_low:i_high), & - i_high - i_low + 1, E) + i_low - 1 - end if - - ! check for rare case where two energy points are the same - if (grid % energy(i_grid) == grid % energy(i_grid + 1)) & - i_grid = i_grid + 1 - - ! calculate interpolation factor - f = (E - grid % energy(i_grid)) / & - (grid % energy(i_grid + 1) - grid % energy(i_grid)) - end associate - - nu_photon_total = ZERO - - ! Calculate nu-photon total cross section - do m = 1, size(this % reactions) - associate (rx => this % reactions(m)) - j = rx % xs(t) % threshold - do k = 1, size(rx % products) - if (rx % products(k) % particle == PHOTON) then - if (i_grid >= j) then - rx_xs = (ONE - f) * rx % xs(t) % value(i_grid - j + 1) & - + f * rx % xs(t) % value(i_grid - j + 2) - nu_photon_total = nu_photon_total + rx_xs * & - rx % products(k) % yield % evaluate(E) - end if - end if - end do - end associate - end do - - end function compute_nuclide_nu_photon_total - subroutine nuclide_init_grid(this, E_min, E_max, M) class(Nuclide), intent(inout) :: this real(8), intent(in) :: E_min ! Minimum energy in MeV diff --git a/src/output.F90 b/src/output.F90 index 479fe6a496..14e69ac510 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -751,12 +751,18 @@ contains score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity" score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power" score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power" - score_names(abs(SCORE_FISS_Q_FRAGMENTS)) = "Fragment fission power" - score_names(abs(SCORE_FISS_Q_BETAS)) = "Beta fission power" - score_names(abs(SCORE_Q_ELASTIC)) = "Elastic scatter power" - score_names(abs(SCORE_Q_PHOTONS)) = "Photon energy deposition" - score_names(abs(SCORE_Q_ELECTRONS)) = "Electron energy deposition" - score_names(abs(SCORE_Q_POSITRONS)) = "Positron energy deposition" + score_names(abs(SCORE_FISS_Q_PROMPT_NEUTRONS)) = "Prompt neutron power" + score_names(abs(SCORE_FISS_Q_DELAYED_NEUTRONS)) = "Delayed neutron power" + score_names(abs(SCORE_FISS_Q_FRAGMENTS)) = "Fission fragment power" + score_names(abs(SCORE_FISS_Q_BETAS)) = "Fission betas power" + score_names(abs(SCORE_FISS_Q_PROMPT_PHOTONS)) = "Prompt photon power" + score_names(abs(SCORE_FISS_Q_DELAYED_PHOTONS)) = "Delayed photon power" + score_names(abs(SCORE_FISS_Q_NEUTRINOS)) = "Fission neutrino power" + score_names(abs(SCORE_Q_PHOTONS)) = "Photon power" + score_names(abs(SCORE_Q_ELECTRONS)) = "Electron power" + score_names(abs(SCORE_Q_POSITRONS)) = "Positron power" + score_names(abs(SCORE_Q_ELASTIC)) = "Elastic scattering power" + score_names(abs(SCORE_HEATING)) = "Heating power" ! Create filename for tally output filename = trim(path_output) // "tallies.out" @@ -910,6 +916,7 @@ contains indent = indent + 2 k = 0 + do l = 1, t % n_user_score_bins k = k + 1 score_index = score_index + 1 @@ -965,6 +972,7 @@ contains end select end associate end do + indent = indent - 2 end do diff --git a/src/physics.F90 b/src/physics.F90 index 96903f54b6..4bc4f609b0 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -112,7 +112,9 @@ contains ! 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 @@ -522,43 +524,50 @@ contains 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 - type(Nuclide), pointer :: nuc ! Get pointer to nuclide - nuc => nuclides(i_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) % nu_photon_total - prob = ZERO + ! 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 + ! 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 + threshold = rx % xs(i_temp) % threshold - ! if energy is below threshold for this reaction, skip it - if (i_grid < threshold) cycle + ! 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 + ! 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) exit REACTION_LOOP - end if - end do - end associate - end do REACTION_LOOP + 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 @@ -1649,7 +1658,6 @@ contains integer :: i_reaction ! index in nuc % reactions array integer :: i_product ! index in nuc % reactions % products array - type(Reaction), pointer :: rx real(8) :: nu_t real(8) :: mu @@ -1659,7 +1667,7 @@ contains integer :: i ! Sample the number of photons produced - nu_t = micro_xs(i_nuclide) % nu_photon_total / micro_xs(i_nuclide) % total + 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 @@ -1671,10 +1679,10 @@ contains ! Sample the reaction and product call sample_photon_product(i_nuclide, p % E, i_reaction, i_product) - rx => nuclides(i_nuclide) % reactions(i_reaction) ! Sample the outgoing energy and angle - call rx % products(i_product) % sample(p % E, E, mu) + 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) diff --git a/src/tally.F90 b/src/tally.F90 index 85bbaea807..43cd7068a8 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,8 +118,6 @@ contains case (SCORE_FLUX, SCORE_FLUX_YN) - if (p % type /= NEUTRON) cycle SCORE_LOOP - 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 @@ -130,7 +129,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 @@ -140,8 +144,6 @@ contains case (SCORE_TOTAL, SCORE_TOTAL_YN) - if (p % type /= NEUTRON) cycle SCORE_LOOP - if (t % estimator == ESTIMATOR_ANALOG) then ! All events will score to the total reaction rate. We can just ! use the weight of the particle entering the collision as the @@ -165,8 +167,6 @@ contains case (SCORE_INVERSE_VELOCITY) - if (p % type /= NEUTRON) cycle SCORE_LOOP - ! make sure the correct energy is used if (t % estimator == ESTIMATOR_TRACKLENGTH) then E = p % E @@ -200,8 +200,6 @@ contains case (SCORE_SCATTER, SCORE_SCATTER_N) - if (p % type /= NEUTRON) cycle SCORE_LOOP - if (t % estimator == ESTIMATOR_ANALOG) then ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP @@ -223,8 +221,6 @@ contains case (SCORE_SCATTER_PN) - if (p % type /= NEUTRON) cycle SCORE_LOOP - ! Only analog estimators are available. ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) then @@ -239,8 +235,6 @@ contains case (SCORE_SCATTER_YN) - if (p % type /= NEUTRON) cycle SCORE_LOOP - ! Only analog estimators are available. ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) then @@ -255,8 +249,6 @@ contains case (SCORE_NU_SCATTER, SCORE_NU_SCATTER_N) - if (p % type /= NEUTRON) cycle SCORE_LOOP - ! Only analog estimators are available. ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP @@ -282,8 +274,6 @@ contains case (SCORE_NU_SCATTER_PN) - if (p % type /= NEUTRON) cycle SCORE_LOOP - ! Only analog estimators are available. ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) then @@ -312,8 +302,6 @@ contains case (SCORE_NU_SCATTER_YN) - if (p % type /= NEUTRON) cycle SCORE_LOOP - ! Only analog estimators are available. ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) then @@ -342,8 +330,6 @@ contains case (SCORE_ABSORPTION) - if (p % type /= NEUTRON) cycle SCORE_LOOP - if (t % estimator == ESTIMATOR_ANALOG) then if (survival_biasing) then ! No absorption events actually occur if survival biasing is on -- @@ -368,7 +354,7 @@ contains case (SCORE_FISSION) - if (p % type /= NEUTRON) cycle SCORE_LOOP + if (material_xs % absorption == ZERO) cycle SCORE_LOOP if (t % estimator == ESTIMATOR_ANALOG) then if (survival_biasing) then @@ -402,7 +388,7 @@ contains case (SCORE_NU_FISSION) - if (p % type /= NEUTRON) cycle SCORE_LOOP + if (material_xs % absorption == ZERO) cycle SCORE_LOOP if (t % estimator == ESTIMATOR_ANALOG) then if (survival_biasing .or. p % fission) then @@ -449,7 +435,7 @@ contains case (SCORE_PROMPT_NU_FISSION) - if (p % type /= NEUTRON) cycle SCORE_LOOP + if (material_xs % absorption == ZERO) cycle SCORE_LOOP ! make sure the correct energy is used if (t % estimator == ESTIMATOR_TRACKLENGTH) then @@ -522,7 +508,7 @@ contains case (SCORE_DELAYED_NU_FISSION) - if (p % type /= NEUTRON) cycle SCORE_LOOP + if (material_xs % absorption == ZERO) cycle SCORE_LOOP ! make sure the correct energy is used if (t % estimator == ESTIMATOR_TRACKLENGTH) then @@ -723,7 +709,7 @@ contains case (SCORE_DECAY_RATE) - if (p % type /= NEUTRON) cycle SCORE_LOOP + if (material_xs % absorption == ZERO) cycle SCORE_LOOP ! make sure the correct energy is used if (t % estimator == ESTIMATOR_TRACKLENGTH) then @@ -1012,7 +998,7 @@ contains case (SCORE_KAPPA_FISSION) - if (p % type /= NEUTRON) cycle SCORE_LOOP + 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 @@ -1082,15 +1068,11 @@ contains case (SCORE_EVENTS) - if (p % type /= NEUTRON) cycle SCORE_LOOP - ! Simply count number of scoring events score = ONE case (ELASTIC) - if (p % type /= NEUTRON) cycle SCORE_LOOP - if (t % estimator == ESTIMATOR_ANALOG) then ! Check if event MT matches if (p % event_MT /= ELASTIC) cycle SCORE_LOOP @@ -1104,9 +1086,12 @@ contains end if end if - case (SCORE_FISS_Q_PROMPT) + case (SCORE_FISS_Q_PROMPT, SCORE_FISS_Q_RECOV, SCORE_FISS_Q_FRAGMENTS, & + SCORE_FISS_Q_PROMPT_NEUTRONS, SCORE_FISS_Q_DELAYED_NEUTRONS, & + SCORE_FISS_Q_PROMPT_PHOTONS, SCORE_FISS_Q_DELAYED_PHOTONS, & + SCORE_FISS_Q_NEUTRINOS, SCORE_FISS_Q_BETAS) - if (p % type /= NEUTRON) cycle SCORE_LOOP + if (material_xs % absorption == ZERO) cycle SCORE_LOOP score = ZERO @@ -1118,10 +1103,29 @@ 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) + else if (score_bin == SCORE_FISS_Q_FRAGMENTS) then + xs = nuc % fission_q_fragments % evaluate(p % last_E) + else if (score_bin == SCORE_FISS_Q_PROMPT_NEUTRONS) then + xs = nuc % fission_q_prompt_neutrons % evaluate(p % last_E) + else if (score_bin == SCORE_FISS_Q_DELAYED_NEUTRONS) then + xs = nuc % fission_q_delayed_neutrons % evaluate(p % last_E) + else if (score_bin == SCORE_FISS_Q_PROMPT_PHOTONS) then + xs = nuc % fission_q_prompt_photons % evaluate(p % last_E) + else if (score_bin == SCORE_FISS_Q_DELAYED_PHOTONS) then + xs = nuc % fission_q_delayed_photons % evaluate(p % last_E) + else if (score_bin == SCORE_FISS_Q_NEUTRINOS) then + xs = nuc % fission_q_neutrinos % evaluate(p % last_E) + else if (score_bin == SCORE_FISS_Q_BETAS) then + xs = nuc % fission_q_betas % 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 @@ -1131,263 +1135,236 @@ contains ! particle's weight entering the collision as the estimate for ! the fission energy production rate associate (nuc => nuclides(p % event_nuclide)) + if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & + 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) + else if (score_bin == SCORE_FISS_Q_FRAGMENTS) then + xs = nuc % fission_q_fragments % evaluate(p % last_E) + else if (score_bin == SCORE_FISS_Q_PROMPT_NEUTRONS) then + xs = nuc % fission_q_prompt_neutrons % evaluate(p % last_E) + else if (score_bin == SCORE_FISS_Q_DELAYED_NEUTRONS) then + xs = nuc % fission_q_delayed_neutrons % evaluate(p % last_E) + else if (score_bin == SCORE_FISS_Q_PROMPT_PHOTONS) then + xs = nuc % fission_q_prompt_photons % evaluate(p % last_E) + else if (score_bin == SCORE_FISS_Q_DELAYED_PHOTONS) then + xs = nuc % fission_q_delayed_photons % evaluate(p % last_E) + else if (score_bin == SCORE_FISS_Q_NEUTRINOS) then + xs = nuc % fission_q_neutrinos % evaluate(p % last_E) + else if (score_bin == SCORE_FISS_Q_BETAS) then + xs = nuc % fission_q_betas % evaluate(p % last_E) + end if + + score = p % last_wgt * xs * flux & + * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption + 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 + associate (nuc => nuclides(i_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 + + 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) + else if (score_bin == SCORE_FISS_Q_FRAGMENTS) then + xs = nuc % fission_q_fragments % evaluate(E) + else if (score_bin == SCORE_FISS_Q_PROMPT_NEUTRONS) then + xs = nuc % fission_q_prompt_neutrons % evaluate(E) + else if (score_bin == SCORE_FISS_Q_DELAYED_NEUTRONS) then + xs = nuc % fission_q_delayed_neutrons % evaluate(E) + else if (score_bin == SCORE_FISS_Q_PROMPT_PHOTONS) then + xs = nuc % fission_q_prompt_photons % evaluate(E) + else if (score_bin == SCORE_FISS_Q_DELAYED_PHOTONS) then + xs = nuc % fission_q_delayed_photons % evaluate(E) + else if (score_bin == SCORE_FISS_Q_NEUTRINOS) then + xs = nuc % fission_q_neutrinos % evaluate(E) + else if (score_bin == SCORE_FISS_Q_BETAS) then + xs = nuc % fission_q_betas % evaluate(E) + end if + + score = micro_xs(i_nuclide) % fission * atom_density * flux * xs 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 + + 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) + else if (score_bin == SCORE_FISS_Q_FRAGMENTS) then + xs = nuc % fission_q_fragments % evaluate(E) + else if (score_bin == SCORE_FISS_Q_PROMPT_NEUTRONS) then + xs = nuc % fission_q_prompt_neutrons % evaluate(E) + else if (score_bin == SCORE_FISS_Q_DELAYED_NEUTRONS) then + xs = nuc % fission_q_delayed_neutrons % evaluate(E) + else if (score_bin == SCORE_FISS_Q_PROMPT_PHOTONS) then + xs = nuc % fission_q_prompt_photons % evaluate(E) + else if (score_bin == SCORE_FISS_Q_DELAYED_PHOTONS) then + xs = nuc % fission_q_delayed_photons % evaluate(E) + else if (score_bin == SCORE_FISS_Q_NEUTRINOS) then + xs = nuc % fission_q_neutrinos % evaluate(E) + else if (score_bin == SCORE_FISS_Q_BETAS) then + xs = nuc % fission_q_betas % 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 end if - case (SCORE_FISS_Q_RECOV) + case (SCORE_Q_ELECTRONS) - if (p % type /= NEUTRON) cycle SCORE_LOOP - - 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) & - * 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 - 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_recov)) then - score = micro_xs(i_nuclide) % fission * atom_density * flux & - * nuclides(i_nuclide) % fission_q_recov % 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_recov)) then - score = score + micro_xs(i_nuc) % fission * atom_density_ & - * flux & - * nuclides(i_nuc) % fission_q_recov % evaluate(E) - end if - end do - end if - end if + ! Electron energy deposition + if (p % type == ELECTRON .and. electron_treatment == ELECTRON_LED) then + score = p % last_wgt * p % last_E end if - case (SCORE_FISS_Q_FRAGMENTS) - if (p % type /= NEUTRON) cycle SCORE_LOOP + case (SCORE_Q_POSITRONS) - 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_fragments)) then - score = p % absorb_wgt & - * nuc % fission_q_fragments % evaluate(p % last_E) & - * 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_fragments)) then - score = p % last_wgt & - * nuc % fission_q_fragments % 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_fragments)) then - score = micro_xs(i_nuclide) % fission * atom_density * flux & - * nuclides(i_nuclide) % fission_q_fragments % 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_fragments)) then - score = score + micro_xs(i_nuc) % fission * atom_density_ & - * flux & - * nuclides(i_nuc) % fission_q_fragments % evaluate(E) - end if - end do - end if - end if + ! Positron energy deposition + if (p % type == POSITRON .and. electron_treatment == ELECTRON_LED) then + score = p % last_wgt * p % last_E end if - case (SCORE_FISS_Q_BETAS) + case (SCORE_Q_PHOTONS) - if (p % type /= NEUTRON) cycle SCORE_LOOP - - 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_betas)) then - score = p % absorb_wgt & - * nuc % fission_q_betas % evaluate(p % last_E) & - * 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_betas)) then - score = p % last_wgt & - * nuc % fission_q_betas % 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_betas)) then - score = micro_xs(i_nuclide) % fission * atom_density * flux & - * nuclides(i_nuclide) % fission_q_betas % 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_betas)) then - score = score + micro_xs(i_nuc) % fission * atom_density_ & - * flux & - * nuclides(i_nuc) % fission_q_betas % evaluate(E) - end if - end do - end if - end if + ! Photon energy deposition + if (p % type == PHOTON .and. p % last_E < energy_cutoff(PHOTON)) then + score = p % last_wgt * p % last_E end if case (SCORE_Q_ELASTIC) - if (p % type /= NEUTRON) cycle SCORE_LOOP - - ! Skip any non-elastic scatter events - if (p % event_MT /= ELASTIC) cycle SCORE_LOOP - - score = p % wgt * (p % last_E - p % E) - - case (SCORE_Q_PHOTONS) - - if (p % type /= PHOTON) cycle SCORE_LOOP - - ! Skip if energy above cutoff - if (p % last_E > energy_cutoff(PHOTON)) cycle SCORE_LOOP - - score = p % wgt * p % last_E - - case (SCORE_Q_ELECTRONS) - - if (p % type /= ELECTRON) cycle SCORE_LOOP - - if (electron_treatment == ELECTRON_LED) then - score = p % wgt * p % last_E + ! Elastic scattering + if (p % event_MT == ELASTIC) then + score = p % last_wgt * (p % last_E - p % E) end if - case (SCORE_Q_POSITRONS) - if (p % type /= POSITRON) cycle SCORE_LOOP + case (SCORE_HEATING) - if (electron_treatment == ELECTRON_LED) then - score = p % wgt * p % last_E + ! Elastic scattering + if (p % event_MT == ELASTIC) then + score = p % last_wgt * (p % last_E - p % E) + + ! Photon energy deposition + else if (p % type == PHOTON) then + if(p % last_E < energy_cutoff(PHOTON)) then + score = p % last_wgt * p % last_E + end if + + ! Electron energy deposition + else if (p % type == ELECTRON) then + if(electron_treatment == ELECTRON_LED) then + score = p % last_wgt * p % last_E + end if + + ! Positron energy deposition + else if (p % type == POSITRON) then + if (electron_treatment == ELECTRON_LED) then + score = p % last_wgt * p % last_E + end if + + ! Fission fragments, betas, and gammas (if photon_transport off) + else + + if (material_xs % absorption == ZERO) cycle SCORE_LOOP + + score = ZERO + + 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)) + score = ZERO + + if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & + allocated(nuc % fission_q_betas)) then + + score = score + p % absorb_wgt & + * nuc % fission_q_fragments % evaluate(p % last_E) & + * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption * flux + + score = score + p % absorb_wgt & + * nuc % fission_q_betas % evaluate(p % last_E) & + * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption * flux + + if (.not. photon_transport) then + score = score + p % absorb_wgt & + * nuc % fission_q_prompt_photons % evaluate(p % last_E) & + * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption * flux + + score = score + p % absorb_wgt & + * nuc % fission_q_delayed_photons % evaluate(p % last_E) & + * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption * flux + end if + end if + end associate + else + ! Skip any non-absorption events + if (p % event /= EVENT_ABSORB) 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 (micro_xs(p % event_nuclide) % absorption > ZERO .and. & + allocated(nuc % fission_q_betas)) then + + score = score + p % last_wgt & + * nuc % fission_q_fragments % evaluate(p % last_E) & + * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption * flux + + score = score + p % last_wgt & + * nuc % fission_q_betas % evaluate(p % last_E) & + * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption * flux + + if (.not. photon_transport) then + score = score + p % last_wgt & + * nuc % fission_q_prompt_photons % evaluate(p % last_E) & + * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption * flux + + score = score + p % last_wgt & + * nuc % fission_q_delayed_photons % evaluate(p % last_E) & + * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption * flux + end if + end if + end associate + end if end if case default @@ -3042,6 +3019,7 @@ contains ! for a previous tally. do j = 1, size(t % filter) i_filt = t % filter(j) + if (.not. filter_matches(i_filt) % bins_present) then call filter_matches(i_filt) % bins % clear() call filter_matches(i_filt) % weights % clear() diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index a34b76fc14..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. !=============================================================================== @@ -638,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 1f56462b77..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