From d1b8efc5aa8157fd883e1d2e726f479069f1f31c Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Sun, 2 Jul 2017 13:21:01 -0400 Subject: [PATCH 01/17] implemented sampling of secondary photons from neutron collisons --- src/cross_section.F90 | 30 +++++++--- src/energy_distribution.F90 | 6 +- src/nuclide_header.F90 | 114 ++++++++++++++++++++++++++++++------ src/physics.F90 | 102 ++++++++++++++++++++++++++++++++ 4 files changed, 225 insertions(+), 27 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 918b06ad83..c16b1ee925 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -53,11 +53,12 @@ contains 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 + material_xs % total = ZERO + material_xs % elastic = ZERO + material_xs % absorption = ZERO + material_xs % fission = ZERO + material_xs % nu_fission = ZERO + material_xs % nu_photon_total = ZERO ! Exit subroutine if material is void if (p % material == MATERIAL_VOID) return @@ -137,6 +138,10 @@ contains ! Add contributions to material macroscopic nu-fission cross section material_xs % nu_fission = material_xs % nu_fission + & atom_density * micro_xs(i_nuclide) % nu_fission + + ! Add contributions to material macroscopic nu-fission cross section + material_xs % nu_photon_total = material_xs % nu_photon_total + & + atom_density * micro_xs(i_nuclide) % nu_photon_total end do end associate @@ -257,8 +262,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) % nu_photon_total = ZERO ! Calculate microscopic nuclide total cross section micro_xs(i_nuclide) % total = (ONE - f) * xs % total(i_grid) & @@ -272,6 +278,10 @@ 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) + if (nuc % fissionable) then ! Calculate microscopic nuclide total cross section micro_xs(i_nuclide) % fission = (ONE - f) * xs % fission(i_grid) & @@ -442,6 +452,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 @@ -563,6 +574,11 @@ 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/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/nuclide_header.F90 b/src/nuclide_header.F90 index 9cc2bc7e20..d6ad5f9f6f 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 + use algorithm, only: sort, find, binary_search use constants use dict_header, only: DictIntInt use endf, only: reaction_name, is_fission, is_disappearance @@ -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 :: nu_photon_total(:) ! photon production end type SumXS type :: Nuclide @@ -99,6 +100,7 @@ module nuclide_header 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 @@ -117,6 +119,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 ! Information for S(a,b) use integer :: index_sab ! index in sab_tables (zero means no table) @@ -138,11 +141,12 @@ module nuclide_header !=============================================================================== type MaterialMacroXS - real(8) :: total ! macroscopic total xs - real(8) :: elastic ! macroscopic elastic scattering xs - real(8) :: absorption ! macroscopic absorption xs - real(8) :: fission ! macroscopic fission xs - real(8) :: nu_fission ! macroscopic production xs + real(8) :: total ! macroscopic total xs + real(8) :: elastic ! macroscopic elastic scattering xs + 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 ! Photon cross sections real(8) :: coherent ! macroscopic coherent xs @@ -481,7 +485,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 +505,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)) 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 end do i_fission = 0 @@ -539,6 +545,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) % nu_photon_total(l+j-1) = & + this % sum_xs(t) % nu_photon_total(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 +583,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 @@ -690,6 +704,72 @@ 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/physics.F90 b/src/physics.F90 index 7582190834..8fc3a501a9 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -113,6 +113,9 @@ contains end if end if + ! Create secondary photons + call sample_secondary_photons(p, i_nuclide) + ! If survival biasing is being used, the following subroutine adjusts the ! weight of the particle. Otherwise, it checks to see if absorption occurs @@ -421,6 +424,59 @@ 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 + real(8) :: f + real(8) :: prob + real(8) :: cutoff + real(8) :: yield + type(Nuclide), pointer :: nuc + + ! Get pointer to nuclide + 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 + + ! 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) exit REACTION_LOOP + end if + end do + end associate + end do REACTION_LOOP + + end subroutine sample_photon_product + !=============================================================================== ! ABSORPTION !=============================================================================== @@ -1498,4 +1554,50 @@ 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 + type(Reaction), pointer :: rx + + 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) % nu_photon_total / 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) + rx => nuclides(i_nuclide) % reactions(i_reaction) + + ! Sample the outgoing energy and angle + call rx % 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 From 20b1113f17b18a55bad59e5461badaf50f6e7f7c Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Sun, 2 Jul 2017 13:29:18 -0400 Subject: [PATCH 02/17] removed macroscopic nu_photon_total --- src/cross_section.F90 | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index c16b1ee925..c3cf4cdf9f 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -58,7 +58,6 @@ contains material_xs % absorption = ZERO material_xs % fission = ZERO material_xs % nu_fission = ZERO - material_xs % nu_photon_total = ZERO ! Exit subroutine if material is void if (p % material == MATERIAL_VOID) return @@ -138,10 +137,6 @@ contains ! Add contributions to material macroscopic nu-fission cross section material_xs % nu_fission = material_xs % nu_fission + & atom_density * micro_xs(i_nuclide) % nu_fission - - ! Add contributions to material macroscopic nu-fission cross section - material_xs % nu_photon_total = material_xs % nu_photon_total + & - atom_density * micro_xs(i_nuclide) % nu_photon_total end do end associate From b4f9badfdb78e6edb8fbf46fbb815fafc91ad4a1 Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Sun, 2 Jul 2017 13:43:39 -0400 Subject: [PATCH 03/17] removed unnecessary adjustment of indent --- src/cross_section.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index c3cf4cdf9f..c0dc3b6852 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -53,11 +53,11 @@ contains 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 + 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 From 5da591d36a442a45b6d2c5f5ddf6931bbaa3f037 Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Mon, 3 Jul 2017 17:01:46 -0400 Subject: [PATCH 04/17] created script for adding photo data to cross_sections.xml and turned on photon transport --- openmc/data/library.py | 23 ++++---- scripts/openmc-get-nndc-data | 9 ++++ scripts/openmc-get-photo-endf71 | 93 +++++++++++++++++++++++++++++++++ src/global.F90 | 2 +- src/photon_physics.F90 | 2 +- src/physics.F90 | 2 +- src/tracking.F90 | 4 +- 7 files changed, 117 insertions(+), 18 deletions(-) create mode 100755 scripts/openmc-get-photo-endf71 diff --git a/openmc/data/library.py b/openmc/data/library.py index 2f33743a6e..5522bdc632 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, it if 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/scripts/openmc-get-nndc-data b/scripts/openmc-get-nndc-data index f1da2241ab..43fa213a0b 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() @@ -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..b2ee3bba2a --- /dev/null +++ b/scripts/openmc-get-photo-endf71 @@ -0,0 +1,93 @@ +#!/usr/bin/env python + +from __future__ import print_function +import os +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 IAEA 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-nds.iaea.org/public/download-endf/ENDF-B-VII.1/' + +# ============================================================================== +# DOWNLOAD FILES FROM IAEA SITE AND GENERATE HDF5 LIBRARY + +# Make photo and ard directories +if not os.path.exists('photo'): + os.mkdir('photo') + +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): + + 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) + + 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) + + 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.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/global.F90 b/src/global.F90 index d9da251bad..6d7b936485 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -120,7 +120,7 @@ module global integer :: n_log_bins ! number of bins for logarithmic grid real(8) :: log_spacing ! spacing on logarithmic grid - logical :: photon_transport = .false. + logical :: photon_transport = .true. ! ============================================================================ ! MULTI-GROUP CROSS SECTION RELATED VARIABLES diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 7b7244fe70..39ae89d2b0 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -327,7 +327,7 @@ contains secondary = elm % shells(i_shell) % transition_subshells(2, i_transition) if (secondary == 0) then - ! Non-radiative trnasition -- Auger/Coster-Kronig effect + ! Non-radiative transition -- Auger/Coster-Kronig effect ! TODO: Create electron ! E_electron = transition_energy(i_transition) diff --git a/src/physics.F90 b/src/physics.F90 index 8fc3a501a9..e6d25cc45b 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1595,7 +1595,7 @@ contains uvw = rotate_angle(p % coord(1) % uvw, mu) ! Create the secondary photon - !call p % create_secondary(uvw, E, PHOTON, run_CE=.true.) + call p % create_secondary(uvw, E, PHOTON, run_CE=.true.) end do end subroutine sample_secondary_photons diff --git a/src/tracking.F90 b/src/tracking.F90 index 7650f598ab..ce513624d8 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -135,7 +135,7 @@ contains ! 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 +166,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 From 10c4113e82f3c943f23f2b16512311eafb929868 Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Mon, 3 Jul 2017 23:00:57 -0400 Subject: [PATCH 05/17] added sampling of electron and positron particle types --- scripts/openmc-get-nndc-data | 2 +- src/constants.F90 | 8 ++- src/global.F90 | 3 +- src/input_xml.F90 | 28 ++++++++++ src/physics.F90 | 104 +++++++++++++++++++++++++++++++---- 5 files changed, 130 insertions(+), 15 deletions(-) diff --git a/scripts/openmc-get-nndc-data b/scripts/openmc-get-nndc-data index 43fa213a0b..8afd98cf70 100755 --- a/scripts/openmc-get-nndc-data +++ b/scripts/openmc-get-nndc-data @@ -113,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 diff --git a/src/constants.F90 b/src/constants.F90 index 53253e23db..f3f73c7ba1 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 :: & @@ -450,6 +451,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/global.F90 b/src/global.F90 index 6d7b936485..680fed2cb9 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 = .true. + 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..26a39de28c 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -267,6 +267,34 @@ 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", 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 + 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) diff --git a/src/physics.F90 b/src/physics.F90 index e6d25cc45b..410f9fcc26 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -46,8 +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) + else if (p % type == ELECTRON) then + call sample_electron_reaction(p) + else if (p % type == POSITRON) then + call sample_positron_reaction(p) end if ! Kill particle if energy falls below cutoff @@ -167,6 +171,8 @@ contains 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 ! Sample element within material i_element = sample_element(p) @@ -219,9 +225,20 @@ 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 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 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. @@ -238,21 +255,84 @@ 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. 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 !=============================================================================== From ae0dd2be143bfffec7e0475b6555e482f7c09645 Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Tue, 4 Jul 2017 17:33:18 -0400 Subject: [PATCH 06/17] added more energy deposition tallies and worked on photoelectric rxn treatment --- src/constants.F90 | 10 +- src/endf.F90 | 12 +++ src/input_xml.F90 | 16 +++ src/nuclide_header.F90 | 30 +++++- src/output.F90 | 6 ++ src/photon_header.F90 | 1 + src/photon_physics.F90 | 60 ++++++----- src/physics.F90 | 17 ++-- src/tally.F90 | 220 ++++++++++++++++++++++++++++++++++++++++- 9 files changed, 333 insertions(+), 39 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index f3f73c7ba1..907a3dbcc8 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 = 24 + integer, parameter :: N_SCORE_TYPES = 30 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -343,7 +343,13 @@ module constants 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_DECAY_RATE = -24 ! delayed neutron precursor decay rate + 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 ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 diff --git a/src/endf.F90 b/src/endf.F90 index f3c6b70899..e34ce36713 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -66,6 +66,18 @@ contains string = "fission-q-prompt" case (SCORE_FISS_Q_RECOV) string = "fission-q-recoverable" + case (SCORE_FISS_Q_FRAGMENTS) + string = "fission-q-fragments" + case (SCORE_FISS_Q_BETAS) + string = "fission-q-betas" + case (SCORE_Q_ELASTIC) + string = "q-elastic" + case (SCORE_Q_PHOTONS) + string = "q-photons" + case (SCORE_Q_ELECTRONS) + string = "q-electrons" + case (SCORE_Q_POSITRONS) + string = "q-positrons" ! Normal ENDF-based reactions case (TOTAL_XS) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 26a39de28c..59359ef115 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3872,6 +3872,22 @@ contains t % score_bins(j) = SCORE_FISS_Q_PROMPT case ('fission-q-recoverable') t % score_bins(j) = SCORE_FISS_Q_RECOV + 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 ('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 ('current') t % score_bins(j) = SCORE_CURRENT t % type = TALLY_SURFACE_CURRENT diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index d6ad5f9f6f..66b190607e 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -92,8 +92,10 @@ 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 + class(Function1D), allocatable :: fission_q_fragments ! fragments + class(Function1D), allocatable :: fission_q_betas ! betas contains procedure :: clear => nuclide_clear @@ -460,6 +462,18 @@ contains fer_dset = open_dataset(fer_group, 'q_recoverable') call this % fission_q_recov % from_hdf5(fer_dset) call close_dataset(fer_dset) + + ! 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) + call close_dataset(fer_dset) + + ! 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) + call close_dataset(fer_dset) else if (temp_str == 'Tabulated1D') then ! Read the prompt Q-value allocate(Tabulated1D :: this % fission_q_prompt) @@ -471,6 +485,18 @@ contains fer_dset = open_dataset(fer_group, 'q_recoverable') call this % fission_q_recov % from_hdf5(fer_dset) call close_dataset(fer_dset) + + ! 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) + call close_dataset(fer_dset) + + ! 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) + call close_dataset(fer_dset) else call fatal_error('Unrecognized fission energy release format.') end if diff --git a/src/output.F90 b/src/output.F90 index 828af717d8..479fe6a496 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -751,6 +751,12 @@ 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" ! Create filename for tally output filename = trim(path_output) // "tallies.out" diff --git a/src/photon_header.F90 b/src/photon_header.F90 index a5cdd9791b..2588f09488 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -209,6 +209,7 @@ contains this % shells(i) % transition_subshells(:,:) = int(matrix(1:2, :), 4) this % shells(i) % transition_energy(:) = matrix(3, :) this % shells(i) % transition_probability(:) = matrix(4, :) + deallocate(matrix) end if call close_dataset(dset_id) diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 39ae89d2b0..3ce47c0da8 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -305,13 +305,22 @@ contains integer :: secondary real(8) :: c real(8) :: rn - real(8) :: E real(8) :: mu real(8) :: phi real(8) :: uvw(3) + real(8) :: E - ! 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 + ! 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) + + ! Get the transition energy + E = elm % shells(i_shell) % transition_energy(i_transition) + + if (secondary /= 0) then ! Non-radiative transition -- Auger/Coster-Kronig effect - ! TODO: Create electron - ! E_electron = transition_energy(i_transition) + ! Create auger electron + call p % create_secondary(uvw, E, ELECTRON, run_ce=.true.) - ! 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 + ! 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 410f9fcc26..88449ca8f6 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -54,13 +54,6 @@ contains call sample_positron_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 - end if - ! Display information about collision if (verbosity >= 10 .or. trace) then if (p % type == NEUTRON) then @@ -174,6 +167,13 @@ contains real(8) :: E_electron ! electron energy real(8) :: uvw(3) ! new direction + ! 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) p % event_nuclide = i_element @@ -242,6 +242,8 @@ contains 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 @@ -267,6 +269,7 @@ contains call p % create_secondary(-uvw, E_electron, POSITRON, .true.) p % event_MT = PAIR_PROD p % alive = .false. + p % E = ZERO end if end subroutine sample_photon_reaction diff --git a/src/tally.F90 b/src/tally.F90 index 14940551b7..85bbaea807 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -115,8 +115,10 @@ contains select case(score_bin) - 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 @@ -137,6 +139,9 @@ 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 @@ -160,6 +165,8 @@ 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 @@ -192,6 +199,9 @@ 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 @@ -212,6 +222,9 @@ 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 @@ -225,6 +238,9 @@ 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 @@ -238,6 +254,9 @@ 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 @@ -262,6 +281,9 @@ 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 @@ -289,6 +311,9 @@ 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 @@ -316,6 +341,9 @@ 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 -- @@ -339,6 +367,9 @@ contains case (SCORE_FISSION) + + if (p % type /= NEUTRON) 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 +401,15 @@ contains case (SCORE_NU_FISSION) + + if (p % type /= NEUTRON) 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 +448,9 @@ contains case (SCORE_PROMPT_NU_FISSION) + + if (p % type /= NEUTRON) cycle SCORE_LOOP + ! make sure the correct energy is used if (t % estimator == ESTIMATOR_TRACKLENGTH) then E = p % E @@ -484,6 +522,8 @@ contains case (SCORE_DELAYED_NU_FISSION) + if (p % type /= NEUTRON) cycle SCORE_LOOP + ! make sure the correct energy is used if (t % estimator == ESTIMATOR_TRACKLENGTH) then E = p % E @@ -683,6 +723,8 @@ contains case (SCORE_DECAY_RATE) + if (p % type /= NEUTRON) cycle SCORE_LOOP + ! make sure the correct energy is used if (t % estimator == ESTIMATOR_TRACKLENGTH) then E = p % E @@ -969,6 +1011,9 @@ contains end if case (SCORE_KAPPA_FISSION) + + if (p % type /= NEUTRON) 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)) @@ -1036,10 +1081,16 @@ contains end if 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 @@ -1054,6 +1105,9 @@ contains end if case (SCORE_FISS_Q_PROMPT) + + if (p % type /= NEUTRON) cycle SCORE_LOOP + score = ZERO if (t % estimator == ESTIMATOR_ANALOG) then @@ -1114,6 +1168,9 @@ contains end if case (SCORE_FISS_Q_RECOV) + + if (p % type /= NEUTRON) cycle SCORE_LOOP + score = ZERO if (t % estimator == ESTIMATOR_ANALOG) then @@ -1173,6 +1230,166 @@ contains end if end if + case (SCORE_FISS_Q_FRAGMENTS) + + 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_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 + end if + + case (SCORE_FISS_Q_BETAS) + + 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 + 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 + end if + + case (SCORE_Q_POSITRONS) + + if (p % type /= POSITRON) cycle SCORE_LOOP + + if (electron_treatment == ELECTRON_LED) then + score = p % wgt * p % last_E + end if + case default if (t % estimator == ESTIMATOR_ANALOG) then ! Any other score is assumed to be a MT number. Thus, we just need @@ -1378,6 +1595,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 From bb61ff78d270b44cb9731af2cec91e14693dcd3a Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Tue, 4 Jul 2017 18:33:33 -0400 Subject: [PATCH 07/17] fixed issue with sampling distance to collison for electrons and positrons --- src/cross_section.F90 | 25 +++++++++++-------------- src/tracking.F90 | 5 +++-- 2 files changed, 14 insertions(+), 16 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index c0dc3b6852..817a874445 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 @@ -595,13 +599,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/tracking.F90 b/src/tracking.F90 index ce513624d8..1f56462b77 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -114,7 +114,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,7 +135,6 @@ contains call score_tracklength_tally(p, distance) end if - ! Score track-length estimate of k-eff if (run_mode == MODE_EIGENVALUE .and. p % type == NEUTRON) then global_tally_tracklength = global_tally_tracklength + p % wgt * & From 58db9d4a08ffb4b51e2e7a0871916a259fed60d0 Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Tue, 4 Jul 2017 21:40:38 -0400 Subject: [PATCH 08/17] turned on samping of secondary photons only if photon transport on --- src/physics.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/physics.F90 b/src/physics.F90 index 88449ca8f6..96903f54b6 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -111,7 +111,9 @@ contains end if ! Create secondary photons - call sample_secondary_photons(p, i_nuclide) + if (photon_transport) then + call sample_secondary_photons(p, i_nuclide) + 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 From 064ff464aff7210e408bd37418c412e9f07e0a1b Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Tue, 4 Jul 2017 23:10:26 -0400 Subject: [PATCH 09/17] fixed issues with using a mesh filter and reading in U234 xs --- src/nuclide_header.F90 | 45 +++++++++++++++++++++++++++++++----------- src/tally_filter.F90 | 6 ++++++ 2 files changed, 39 insertions(+), 12 deletions(-) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 66b190607e..716c482067 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -452,6 +452,7 @@ contains 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) @@ -463,18 +464,8 @@ contains call this % fission_q_recov % from_hdf5(fer_dset) call close_dataset(fer_dset) - ! 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) - call close_dataset(fer_dset) - - ! 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) - 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) @@ -485,12 +476,42 @@ contains 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.') + end if + + fer_dset = open_dataset(fer_group, 'fragments') + 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) + 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) call close_dataset(fer_dset) + else + call fatal_error('Unrecognized fission fragment energy release format.') + end if + + fer_dset = open_dataset(fer_group, 'betas') + 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) + call close_dataset(fer_dset) + + else if (temp_str == 'Tabulated1D') then ! Read the beta energy Q-value allocate(Tabulated1D :: this % fission_q_betas) @@ -498,7 +519,7 @@ contains call this % fission_q_betas % from_hdf5(fer_dset) call close_dataset(fer_dset) else - call fatal_error('Unrecognized fission energy release format.') + call fatal_error('Unrecognized beta energy release format.') end if call close_group(fer_group) end if diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index dfb89d4330..a34b76fc14 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -313,6 +313,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. From 836a9836098b0fc79097ae96805b40abf6fc2063 Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Wed, 5 Jul 2017 09:05:23 -0400 Subject: [PATCH 10/17] fixed typo in library.py --- openmc/data/library.py | 2 +- src/photon_header.F90 | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/openmc/data/library.py b/openmc/data/library.py index 5522bdc632..baf8f6ae0b 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -71,7 +71,7 @@ class DataLibrary(EqualityMixin): path : str Path to file to write. Defaults to 'cross_sections.xml'. append : bool - Whether to append to an existing file, it if exists. + Whether to append to an existing file, if it exists. Defaults to False. """ diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 2588f09488..a5cdd9791b 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -209,7 +209,6 @@ contains this % shells(i) % transition_subshells(:,:) = int(matrix(1:2, :), 4) this % shells(i) % transition_energy(:) = matrix(3, :) this % shells(i) % transition_probability(:) = matrix(4, :) - deallocate(matrix) end if call close_dataset(dset_id) From 8da0636ab827b1d409425080aa81ef215b738be6 Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Thu, 13 Jul 2017 09:32:02 -0400 Subject: [PATCH 11/17] 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 From dc32bccfc35b476045d6f158c731e186c7fc0538 Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Thu, 13 Jul 2017 09:39:26 -0400 Subject: [PATCH 12/17] fixed error in issuing particle filter warning --- src/input_xml.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 34b14cddec..82cbbe6632 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -4170,10 +4170,10 @@ contains 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 + 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 % bins % data(l))) // " will have no scores") + &type " // trim(to_str(filt % particles(l))) // " will have no scores") end if end do end select From dc71fa291538507d35261c17abb07ff176c90cd7 Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Thu, 13 Jul 2017 09:55:55 -0400 Subject: [PATCH 13/17] fixed assorted issues with heating scores --- scripts/openmc-get-photo-endf71 | 2 +- src/constants.F90 | 22 +++++++++++----------- src/cross_section.F90 | 6 +++--- src/endf.F90 | 6 ++++++ src/global.F90 | 1 - src/input_xml.F90 | 2 +- src/output.F90 | 2 -- src/photon_physics.F90 | 2 +- 8 files changed, 23 insertions(+), 20 deletions(-) diff --git a/scripts/openmc-get-photo-endf71 b/scripts/openmc-get-photo-endf71 index 7fd06a0a24..9920b27e28 100755 --- a/scripts/openmc-get-photo-endf71 +++ b/scripts/openmc-get-photo-endf71 @@ -14,7 +14,7 @@ import openmc.data from openmc.data import ATOMIC_SYMBOL description = """ -Download ENDF/B-VII.1 ENDF data from IAEA for photo-atomic and atomic +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. diff --git a/src/constants.F90 b/src/constants.F90 index c04ecf1b2b..bb06497229 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -343,18 +343,18 @@ module constants SCORE_INVERSE_VELOCITY = -21, & ! flux-weighted inverse velocity SCORE_HEATING = -22, & ! prompt fission Q-value SCORE_FISS_Q_RECOV = -23, & ! recoverable fission Q-value - SCORE_FISS_Q_PROMPT = -24, & ! recoverable fission Q-value - SCORE_FISS_Q_PROMPT_NEUTRONS = -25, & ! recoverable fission Q-value + SCORE_FISS_Q_PROMPT = -24, & ! prompt 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_FISS_Q_FRAGMENTS = -27, & ! fission fragments Q-value + SCORE_FISS_Q_BETAS = -28, & ! fission betas Q-value + SCORE_FISS_Q_PROMPT_PHOTONS = -29, & ! fission prompt photons Q-value + SCORE_FISS_Q_DELAYED_PHOTONS = -30, & ! fission delayed phtons Q-value + SCORE_FISS_Q_NEUTRINOS = -31, & ! fission neutrinos Q-value + SCORE_Q_PHOTONS = -32, & ! photon energy deposition + SCORE_Q_ELECTRONS = -33, & ! electron energy deposition + SCORE_Q_POSITRONS = -34, & ! positron energy deposition + SCORE_Q_ELASTIC = -35, & ! elastic scattering energy deposition SCORE_DECAY_RATE = -36 ! delayed neutron precursor decay rate ! Maximum scattering order supported diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 325f1ed841..37a3768455 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -261,9 +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) % photon_prod = 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) & diff --git a/src/endf.F90 b/src/endf.F90 index e34ce36713..f67b231c11 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -70,6 +70,12 @@ contains string = "fission-q-fragments" case (SCORE_FISS_Q_BETAS) string = "fission-q-betas" + case (SCORE_FISS_Q_PROMPT_PHOTONS) + string = "fission-q-prompt-photons" + case (SCORE_FISS_Q_DELAYED_PHOTONS) + string = "fission-q-delayed-photons" + case (SCORE_FISS_Q_NEUTRINOS) + string = "fission-q-neutrinos" case (SCORE_Q_ELASTIC) string = "q-elastic" case (SCORE_Q_PHOTONS) diff --git a/src/global.F90 b/src/global.F90 index 28ff57cc60..30d824a6d2 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -59,7 +59,6 @@ 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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 82cbbe6632..6f3c9be71a 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -285,7 +285,7 @@ contains if (check_for_node(root, "photon_transport")) then call get_node_value(root, "photon_transport", photon_transport) - if (.not. run_CE) then + if (.not. run_CE .and. photon_transport) then call fatal_error("Photon transport is not currently supported & &in Multi-group mode") end if diff --git a/src/output.F90 b/src/output.F90 index 14e69ac510..7acd303fda 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -916,7 +916,6 @@ contains indent = indent + 2 k = 0 - do l = 1, t % n_user_score_bins k = k + 1 score_index = score_index + 1 @@ -972,7 +971,6 @@ contains end select end associate end do - indent = indent - 2 end do diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 3ce47c0da8..7d4809803d 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -305,10 +305,10 @@ contains integer :: secondary real(8) :: c real(8) :: rn + real(8) :: E real(8) :: mu real(8) :: phi real(8) :: uvw(3) - real(8) :: E ! If no transitions, assume fluorescent photon from captured free electron if (elm % shells(i_shell) % n_transitions == 0) then From ddda797c0612f9d449971ff799bbde3d7d7664df Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Thu, 13 Jul 2017 10:06:42 -0400 Subject: [PATCH 14/17] removed unnecessary PR edits --- src/constants.F90 | 4 ++-- src/nuclide_header.F90 | 22 +++++++++++----------- src/tally.F90 | 16 +--------------- 3 files changed, 14 insertions(+), 28 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index bb06497229..fa792717b0 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -344,8 +344,8 @@ module constants SCORE_HEATING = -22, & ! prompt fission Q-value SCORE_FISS_Q_RECOV = -23, & ! recoverable fission Q-value SCORE_FISS_Q_PROMPT = -24, & ! prompt 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_PROMPT_NEUTRONS = -25, & ! fission prompt neutrons Q-value + SCORE_FISS_Q_DELAYED_NEUTRONS = -26, & ! fission delayed neutrons Q-value SCORE_FISS_Q_FRAGMENTS = -27, & ! fission fragments Q-value SCORE_FISS_Q_BETAS = -28, & ! fission betas Q-value SCORE_FISS_Q_PROMPT_PHOTONS = -29, & ! fission prompt photons Q-value diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 84de722a7d..276517c4fb 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -96,11 +96,11 @@ 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 + class(Function1D), allocatable :: fission_q_neutrinos ! neutrinos + class(Function1D), allocatable :: fission_q_delayed_neutrons ! delayed neutrons + class(Function1D), allocatable :: fission_q_prompt_neutrons ! prompt neutrons + class(Function1D), allocatable :: fission_q_delayed_photons ! delayed photons + class(Function1D), allocatable :: fission_q_prompt_photons ! prompt photons contains procedure :: clear => nuclide_clear @@ -147,12 +147,12 @@ module nuclide_header !=============================================================================== type MaterialMacroXS - real(8) :: total ! macroscopic total xs - real(8) :: elastic ! macroscopic elastic scattering xs - 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 + real(8) :: total ! macroscopic total xs + real(8) :: elastic ! macroscopic elastic scattering xs + 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 diff --git a/src/tally.F90 b/src/tally.F90 index 43cd7068a8..e26e18f169 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -116,6 +116,7 @@ contains select case(score_bin) + case (SCORE_FLUX, SCORE_FLUX_YN) if (t % estimator == ESTIMATOR_ANALOG) then @@ -143,7 +144,6 @@ contains case (SCORE_TOTAL, SCORE_TOTAL_YN) - 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 @@ -166,7 +166,6 @@ contains case (SCORE_INVERSE_VELOCITY) - ! make sure the correct energy is used if (t % estimator == ESTIMATOR_TRACKLENGTH) then E = p % E @@ -199,7 +198,6 @@ contains case (SCORE_SCATTER, SCORE_SCATTER_N) - if (t % estimator == ESTIMATOR_ANALOG) then ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP @@ -220,7 +218,6 @@ contains case (SCORE_SCATTER_PN) - ! Only analog estimators are available. ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) then @@ -234,7 +231,6 @@ contains case (SCORE_SCATTER_YN) - ! Only analog estimators are available. ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) then @@ -248,7 +244,6 @@ contains case (SCORE_NU_SCATTER, SCORE_NU_SCATTER_N) - ! Only analog estimators are available. ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP @@ -273,7 +268,6 @@ contains case (SCORE_NU_SCATTER_PN) - ! Only analog estimators are available. ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) then @@ -301,7 +295,6 @@ contains case (SCORE_NU_SCATTER_YN) - ! Only analog estimators are available. ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) then @@ -329,7 +322,6 @@ contains case (SCORE_ABSORPTION) - if (t % estimator == ESTIMATOR_ANALOG) then if (survival_biasing) then ! No absorption events actually occur if survival biasing is on -- @@ -1067,12 +1059,10 @@ contains end if case (SCORE_EVENTS) - ! Simply count number of scoring events score = ONE case (ELASTIC) - if (t % estimator == ESTIMATOR_ANALOG) then ! Check if event MT matches if (p % event_MT /= ELASTIC) cycle SCORE_LOOP @@ -1244,7 +1234,6 @@ contains score = p % last_wgt * p % last_E end if - case (SCORE_Q_POSITRONS) ! Positron energy deposition @@ -1266,7 +1255,6 @@ contains score = p % last_wgt * (p % last_E - p % E) end if - case (SCORE_HEATING) ! Elastic scattering @@ -1572,7 +1560,6 @@ 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 @@ -3019,7 +3006,6 @@ 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() From ba006db40f4c75ad9674f3f25c774c628a9c65ba Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Thu, 13 Jul 2017 14:37:51 -0400 Subject: [PATCH 15/17] added non-relativistic Sauter distribution for photoelectron outgoing angle --- src/physics.F90 | 51 +++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 49 insertions(+), 2 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 4bc4f609b0..3848a79a28 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -163,6 +163,7 @@ 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 @@ -170,6 +171,7 @@ contains 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 @@ -231,8 +233,20 @@ contains if (prob > cutoff) then E_electron = p % E - elm % shells(i_shell) % binding_energy - ! Sample angle isotropically - mu = TWO*prn() - ONE + ! 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_MU1: 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_MU1 + end if + end do SAMPLE_MU1 + phi = TWO*PI*prn() uvw(1) = mu uvw(2) = sqrt(ONE - mu*mu)*cos(phi) @@ -251,6 +265,39 @@ contains return end if end do + + ! If no shell was sampled, give the whole phton energy to the electron. + ! See Eqn 3.9 in "Implementing a photon physics model in Serpent 2" by + ! Toni Kaltiaisenaho + + ! 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_MU2: 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_MU2 + end if + end do SAMPLE_MU2 + + 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, p % E, ELECTRON, run_CE=.true.) + + ! TODO: figure out correct event_MT + p % event_MT = 533 + + p % alive = .false. + p % E = ZERO + return end if prob = prob_after end associate From 41272a3aa4f616ee241f5adce4afb10bf5c200b3 Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Thu, 13 Jul 2017 15:58:08 -0400 Subject: [PATCH 16/17] removed option for giving photoelectron entire photon energy --- src/physics.F90 | 44 ++++++++------------------------------------ 1 file changed, 8 insertions(+), 36 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 3848a79a28..9a07764152 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -206,6 +206,11 @@ contains prob = prob + micro_photon_xs(i_element) % incoherent if (prob > cutoff) then call compton_scatter(elm, alpha, alpha_out, mu, .true.) + + ! Create secondary electron + + + p % E = alpha_out*MASS_ELECTRON p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) p % event_MT = INCOHERENT @@ -236,16 +241,16 @@ contains ! 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_MU1: do + 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_MU1 + exit SAMPLE_MU end if - end do SAMPLE_MU1 + end do SAMPLE_MU phi = TWO*PI*prn() uvw(1) = mu @@ -265,39 +270,6 @@ contains return end if end do - - ! If no shell was sampled, give the whole phton energy to the electron. - ! See Eqn 3.9 in "Implementing a photon physics model in Serpent 2" by - ! Toni Kaltiaisenaho - - ! 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_MU2: 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_MU2 - end if - end do SAMPLE_MU2 - - 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, p % E, ELECTRON, run_CE=.true.) - - ! TODO: figure out correct event_MT - p % event_MT = 533 - - p % alive = .false. - p % E = ZERO - return end if prob = prob_after end associate From beb2c40a060b326e5cb566ccdf40f0c4c80b1455 Mon Sep 17 00:00:00 2001 From: samuelshaner Date: Thu, 13 Jul 2017 16:06:05 -0400 Subject: [PATCH 17/17] removed new energy deposition scores --- src/constants.F90 | 20 +---- src/endf.F90 | 18 ---- src/input_xml.F90 | 26 ------ src/nuclide_header.F90 | 116 +------------------------ src/output.F90 | 12 --- src/physics.F90 | 5 -- src/tally.F90 | 189 +---------------------------------------- 7 files changed, 7 insertions(+), 379 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index fa792717b0..16d266a5a1 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 = 36 + integer, parameter :: N_SCORE_TYPES = 24 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -341,21 +341,9 @@ 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_HEATING = -22, & ! prompt fission Q-value - SCORE_FISS_Q_RECOV = -23, & ! recoverable fission Q-value - SCORE_FISS_Q_PROMPT = -24, & ! prompt fission Q-value - SCORE_FISS_Q_PROMPT_NEUTRONS = -25, & ! fission prompt neutrons Q-value - SCORE_FISS_Q_DELAYED_NEUTRONS = -26, & ! fission delayed neutrons Q-value - SCORE_FISS_Q_FRAGMENTS = -27, & ! fission fragments Q-value - SCORE_FISS_Q_BETAS = -28, & ! fission betas Q-value - SCORE_FISS_Q_PROMPT_PHOTONS = -29, & ! fission prompt photons Q-value - SCORE_FISS_Q_DELAYED_PHOTONS = -30, & ! fission delayed phtons Q-value - SCORE_FISS_Q_NEUTRINOS = -31, & ! fission neutrinos Q-value - SCORE_Q_PHOTONS = -32, & ! photon energy deposition - SCORE_Q_ELECTRONS = -33, & ! electron energy deposition - SCORE_Q_POSITRONS = -34, & ! positron energy deposition - SCORE_Q_ELASTIC = -35, & ! elastic scattering energy deposition - SCORE_DECAY_RATE = -36 ! delayed neutron precursor decay rate + 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 integer, parameter :: MAX_ANG_ORDER = 10 diff --git a/src/endf.F90 b/src/endf.F90 index f67b231c11..f3c6b70899 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -66,24 +66,6 @@ contains string = "fission-q-prompt" case (SCORE_FISS_Q_RECOV) string = "fission-q-recoverable" - case (SCORE_FISS_Q_FRAGMENTS) - string = "fission-q-fragments" - case (SCORE_FISS_Q_BETAS) - string = "fission-q-betas" - case (SCORE_FISS_Q_PROMPT_PHOTONS) - string = "fission-q-prompt-photons" - case (SCORE_FISS_Q_DELAYED_PHOTONS) - string = "fission-q-delayed-photons" - case (SCORE_FISS_Q_NEUTRINOS) - string = "fission-q-neutrinos" - case (SCORE_Q_ELASTIC) - string = "q-elastic" - case (SCORE_Q_PHOTONS) - string = "q-photons" - case (SCORE_Q_ELECTRONS) - string = "q-electrons" - case (SCORE_Q_POSITRONS) - string = "q-positrons" ! Normal ENDF-based reactions case (TOTAL_XS) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 6f3c9be71a..5273a21807 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3883,32 +3883,6 @@ 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 ('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 diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 276517c4fb..ab937809ee 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -92,15 +92,8 @@ module nuclide_header ! array; used at tally-time ! Fission energy release - class(Function1D), allocatable :: fission_q_prompt ! fragments and prompt neutrons, gammas - 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 ! neutrinos - class(Function1D), allocatable :: fission_q_delayed_neutrons ! delayed neutrons - class(Function1D), allocatable :: fission_q_prompt_neutrons ! prompt neutrons - class(Function1D), allocatable :: fission_q_delayed_photons ! delayed photons - class(Function1D), allocatable :: fission_q_prompt_photons ! prompt photons + 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 @@ -482,111 +475,6 @@ contains 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 - - ! Q-DELAYED-PHOTONS - fer_dset = open_dataset(fer_group, 'delayed_photons') - call read_attribute(temp_str, fer_dset, 'type') - if (temp_str == 'Polynomial') then - 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 - 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 delayed photon energy release format.') - end if - - ! Q-PROMPT-PHOTONS - fer_dset = open_dataset(fer_group, 'prompt_photons') - call read_attribute(temp_str, fer_dset, 'type') - if (temp_str == 'Polynomial') then - 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 - 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 fission prompt photon energy release format.') - end if - call close_group(fer_group) end if diff --git a/src/output.F90 b/src/output.F90 index 7acd303fda..828af717d8 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -751,18 +751,6 @@ 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_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" diff --git a/src/physics.F90 b/src/physics.F90 index 9a07764152..6923a5386a 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -206,11 +206,6 @@ contains prob = prob + micro_photon_xs(i_element) % incoherent if (prob > cutoff) then call compton_scatter(elm, alpha, alpha_out, mu, .true.) - - ! Create secondary electron - - - p % E = alpha_out*MASS_ELECTRON p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) p % event_MT = INCOHERENT diff --git a/src/tally.F90 b/src/tally.F90 index e26e18f169..3d7ee2cb80 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1076,10 +1076,7 @@ contains end if end if - 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) + case (SCORE_FISS_Q_PROMPT, SCORE_FISS_Q_RECOV) if (material_xs % absorption == ZERO) cycle SCORE_LOOP @@ -1097,20 +1094,6 @@ contains 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 & @@ -1132,20 +1115,6 @@ contains 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 & @@ -1170,20 +1139,6 @@ contains 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 @@ -1202,20 +1157,6 @@ contains 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_ & @@ -1227,134 +1168,6 @@ contains end if end if - case (SCORE_Q_ELECTRONS) - - ! Electron energy deposition - if (p % type == ELECTRON .and. electron_treatment == ELECTRON_LED) then - score = p % last_wgt * p % last_E - end if - - case (SCORE_Q_POSITRONS) - - ! Positron energy deposition - if (p % type == POSITRON .and. electron_treatment == ELECTRON_LED) then - score = p % last_wgt * p % last_E - end if - - case (SCORE_Q_PHOTONS) - - ! 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) - - ! Elastic scattering - if (p % event_MT == ELASTIC) then - score = p % last_wgt * (p % last_E - p % E) - end if - - case (SCORE_HEATING) - - ! 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 if (t % estimator == ESTIMATOR_ANALOG) then ! Any other score is assumed to be a MT number. Thus, we just need