crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

added particle filter and addressed PR comments

Browse source
This commit is contained in:
samuelshaner 2017-07-13 09:32:02 -04:00
parent 836a983609
commit 8da0636ab8
13 changed files with 687 additions and 502 deletions
Download patch file Download diff file Expand all files Collapse all files

40
openmc/filter.py
View file

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

26
openmc/settings.py
View file

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

58
scripts/openmc-get-photo-endf71
View file

@ -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)

32
src/constants.F90
View file

@ -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

13
src/cross_section.F90
View file

@ -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) * &

3
src/global.F90
View file

@ -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
! ============================================================================

118
src/input_xml.F90
View file

@ -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 <scores> specified on tally " &
// trim(to_str(t % id)) // ".")

232
src/nuclide_header.F90
View file

@ -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

20
src/output.F90
View file

@ -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

66
src/physics.F90
View file

@ -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)

506
src/tally.F90
View file

@ -99,6 +99,7 @@ contains
real(8) :: f ! interpolation factor
real(8) :: score ! analog tally score
real(8) :: E ! particle energy
real(8) :: xs ! cross section
i = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
@ -117,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()

65
src/tally_filter.F90
View file

@ -58,6 +58,19 @@ module tally_filter
procedure :: initialize => initialize_material
end type MaterialFilter
!===============================================================================
! PARTICLE specifies which particle tally events reside in.
!===============================================================================
type, extends(TallyFilter) :: ParticleFilter
integer, allocatable :: particles(:)
type(DictIntInt) :: map
contains
procedure :: get_next_bin => get_next_bin_particle
procedure :: to_statepoint => to_statepoint_particle
procedure :: text_label => text_label_particle
procedure :: initialize => initialize_particle
end type ParticleFilter
!===============================================================================
! CELLFILTER specifies which geometric cells tally events reside in.
!===============================================================================
@ -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
!===============================================================================

10
src/tracking.F90
View file

@ -11,7 +11,7 @@ module tracking
use particle_header, only: LocalCoord, Particle
use physics, only: collision
use physics_mg, only: collision_mg
use random_lcg, only: prn
use random_lcg, only: prn, prn_set_stream
use string, only: to_str
use tally, only: score_analog_tally, score_tracklength_tally, &
score_collision_tally, score_surface_current, &
@ -69,6 +69,14 @@ contains
if (active_tallies % size() > 0) call zero_flux_derivs()
EVENT_LOOP: do
! Set the random number stream
if (p % type == NEUTRON) then
call prn_set_stream(STREAM_TRACKING)
else
call prn_set_stream(STREAM_PHOTON)
end if
! If the cell hasn't been determined based on the particle's location,
! initiate a search for the current cell. This generally happens at the
! beginning of the history and again for any secondary particles