diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 11214eed2a..615908dc46 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -29,6 +29,11 @@ fileio.o: constants.o fileio.o: global.o fileio.o: string.o +fission.o: constants.o +fission.o: cross_section_header.o +fission.o: error.o +fission.o: search.o + geometry.o: constants.o geometry.o: datatypes.o geometry.o: error.o @@ -113,6 +118,7 @@ physics.o: constants.o physics.o: cross_section_header.o physics.o: endf.o physics.o: error.o +physics.o: fission.o physics.o: geometry.o physics.o: geometry_header.o physics.o: global.o diff --git a/src/OBJECTS b/src/OBJECTS index cdf0fab349..8f15589b5e 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -9,6 +9,7 @@ endf.o \ energy_grid.o \ error.o \ fileio.o \ +fission.o \ geometry.o \ geometry_header.o \ global.o \ diff --git a/src/constants.f90 b/src/constants.f90 index 200a4b71e0..08a72922b6 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -2,6 +2,11 @@ module constants implicit none + ! Versioning numbers + integer, parameter :: VERSION_MAJOR = 0 + integer, parameter :: VERSION_MINOR = 3 + integer, parameter :: VERSION_RELEASE = 2 + ! Physical constants real(8), parameter :: & PI = 3.1415926535898_8, & ! pi @@ -98,53 +103,53 @@ module constants ! Reaction types integer, parameter :: & - TOTAL_XS = 1, & - ELASTIC = 2, & - N_LEVEL = 4, & - MISC = 5, & - N_2ND = 11, & - N_2N = 16, & - N_3N = 17, & - FISSION = 18, & - N_F = 19, & - N_NF = 20, & - N_2NF = 21, & - N_NA = 22, & - N_N3A = 23, & - N_2NA = 24, & - N_3NA = 25, & - N_NP = 28, & - N_N2A = 29, & - N_2N2A = 30, & - N_ND = 32, & - N_NT = 33, & - N_N3HE = 34, & - N_ND2A = 35, & - N_NT2A = 36, & - N_4N = 37, & - N_3NF = 38, & - N_2NP = 41, & - N_3NP = 42, & - N_N2P = 44, & - N_NPA = 45, & - N_N1 = 51, & - N_N40 = 90, & - N_NC = 91, & - N_GAMMA = 102, & - N_P = 103, & - N_D = 104, & - N_T = 105, & - N_3HE = 106, & - N_A = 107, & - N_2A = 108, & - N_3A = 109, & - N_2P = 111, & - N_PA = 112, & - N_T2A = 113, & - N_D2A = 114, & - N_PD = 115, & - N_PT = 116, & - N_DA = 117 + TOTAL_XS = 1, & + ELASTIC = 2, & + N_LEVEL = 4, & + MISC = 5, & + N_2ND = 11, & + N_2N = 16, & + N_3N = 17, & + N_FISSION = 18, & + N_F = 19, & + N_NF = 20, & + N_2NF = 21, & + N_NA = 22, & + N_N3A = 23, & + N_2NA = 24, & + N_3NA = 25, & + N_NP = 28, & + N_N2A = 29, & + N_2N2A = 30, & + N_ND = 32, & + N_NT = 33, & + N_N3HE = 34, & + N_ND2A = 35, & + N_NT2A = 36, & + N_4N = 37, & + N_3NF = 38, & + N_2NP = 41, & + N_3NP = 42, & + N_N2P = 44, & + N_NPA = 45, & + N_N1 = 51, & + N_N40 = 90, & + N_NC = 91, & + N_GAMMA = 102, & + N_P = 103, & + N_D = 104, & + N_T = 105, & + N_3HE = 106, & + N_A = 107, & + N_2A = 108, & + N_3A = 109, & + N_2P = 111, & + N_PA = 112, & + N_T2A = 113, & + N_D2A = 114, & + N_PD = 115, & + N_PT = 116, & + N_DA = 117 ! Tally distinguishers integer, parameter :: & @@ -159,11 +164,6 @@ module constants NU_POLYNOMIAL = 1, & ! Nu values given by polynomial NU_TABULAR = 2 ! Nu values given by tabular distribution - ! Versioning numbers - integer, parameter :: VERSION_MAJOR = 0 - integer, parameter :: VERSION_MINOR = 3 - integer, parameter :: VERSION_RELEASE = 1 - ! Maximum number of words in a single line, length of line, and length of ! single word integer, parameter :: MAX_WORDS = 500 diff --git a/src/cross_section.f90 b/src/cross_section.f90 index 33ec9fa9c4..5aa10ec81c 100644 --- a/src/cross_section.f90 +++ b/src/cross_section.f90 @@ -84,6 +84,9 @@ contains allocate(nuclides(n_nuclides_total)) allocate(sab_tables(n_sab_tables)) + ! allocate array for microscopic cross section cache + allocate(micro_xs(n_nuclides_total)) + ! loop over all nuclides in xsdata call dict_create(temp_dict) @@ -241,25 +244,36 @@ contains ! determine number of energy points NE = NXS(3) - nuc%n_grid = NE + nuc % n_grid = NE - ! allocate storage for arrays - allocate(nuc%energy(NE)) - allocate(nuc%sigma_t(NE)) - allocate(nuc%sigma_a(NE)) - allocate(nuc%sigma_el(NE)) - allocate(nuc%heating(NE)) + ! allocate storage for energy grid and cross section arrays + allocate(nuc % energy(NE)) + allocate(nuc % total(NE)) + allocate(nuc % elastic(NE)) + allocate(nuc % fission(NE)) + allocate(nuc % absorption(NE)) + allocate(nuc % heating(NE)) + + ! initialize cross sections + nuc % total = ZERO + nuc % elastic = ZERO + nuc % fission = ZERO + nuc % absorption = ZERO + nuc % heating = ZERO + + ! Read data from XSS -- only the energy grid, elastic scattering and heating + ! cross section values are actually read from here. The total and absorption + ! cross sections are reconstructed from the partial reaction data. - ! read data from XSS -- right now the total, absorption and elastic - ! scattering are read in to these special arrays, but in reality, it should - ! be necessary to only store elastic scattering and possibly total - ! cross-section for total material xs generation. XSS_index = 1 - nuc%energy = get_real(NE) - nuc%sigma_t = get_real(NE) - nuc%sigma_a = get_real(NE) - nuc%sigma_el = get_real(NE) - nuc%heating = get_real(NE) + nuc % energy = get_real(NE) + + ! Skip total and absorption + XSS_index = XSS_index + 2*NE + + ! Continue reading elastic scattering and heating + nuc % elastic = get_real(NE) + nuc % heating = get_real(NE) end subroutine read_esz @@ -515,6 +529,7 @@ contains integer :: JXS7 ! index of reactions cross-sections in XSS integer :: LXS ! location of cross-section locators integer :: LOCA ! location of cross-section for given MT + integer :: IE ! reaction's starting index on energy grid integer :: NE ! number of energies for reaction type(Reaction), pointer :: rxn => null() @@ -527,34 +542,60 @@ contains ! allocate array of reactions. Add one since we need to include an elastic ! scattering channel - nuc%n_reaction = NMT + 1 - allocate(nuc%reactions(NMT+1)) + nuc % n_reaction = NMT + 1 + allocate(nuc % reactions(NMT+1)) ! Store elastic scattering cross-section on reaction one - rxn => nuc%reactions(1) - rxn%MT = 2 - rxn%Q_value = ZERO - rxn%TY = 1 - rxn%IE = 1 - allocate(rxn%sigma(nuc%n_grid)) - rxn%sigma = nuc%sigma_el + rxn => nuc % reactions(1) + rxn % MT = 2 + rxn % Q_value = ZERO + rxn % TY = 1 + rxn % IE = 1 + allocate(rxn % sigma(nuc % n_grid)) + rxn % sigma = nuc % elastic + + ! Add contribution of elastic scattering to total cross section + nuc % total = nuc % total + nuc % elastic do i = 1, NMT - rxn => nuc%reactions(i+1) + rxn => nuc % reactions(i+1) ! read MT number, Q-value, and neutrons produced - rxn % MT = XSS(LMT+i-1) - rxn % Q_value = XSS(JXS4+i-1) - rxn % TY = XSS(JXS5+i-1) + rxn % MT = XSS(LMT + i - 1) + rxn % Q_value = XSS(JXS4 + i - 1) + rxn % TY = XSS(JXS5 + i - 1) - ! read cross section values - LOCA = XSS(LXS+i-1) - rxn % IE = XSS(JXS7 + LOCA - 1) + ! read starting energy index + LOCA = XSS(LXS + i - 1) + IE = XSS(JXS7 + LOCA - 1) + rxn % IE = IE + + ! read number of energies cross section values NE = XSS(JXS7 + LOCA) - allocate(rxn%sigma(NE)) + allocate(rxn % sigma(NE)) XSS_index = JXS7 + LOCA + 1 rxn % sigma = get_real(NE) + ! Skip redundant reactions -- this includes total inelastic level + ! scattering, gas production cross sections (MT=200+), and (n,p), (n,d), + ! etc. reactions leaving the nucleus in an excited state + if (rxn % MT == N_LEVEL .or. rxn % MT > N_DA) cycle + + ! Add contribution to total cross section + nuc % total(IE:IE+NE-1) = nuc % total(IE:IE+NE-1) + rxn % sigma + + ! Add contribution to absorption cross section + if (rxn % MT >= N_GAMMA .and. rxn % MT <= N_DA) then + nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma + end if + + ! Add contribution to fission cross section + if (rxn % MT == N_FISSION .or. rxn % MT == N_F .or. rxn % MT == N_NF & + .or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) then + nuc % fissionable = .true. + nuc % fission(IE:IE+NE-1) = nuc % fission(IE:IE+NE-1) + rxn % sigma + end if + ! set defaults rxn % has_angle_dist = .false. rxn % has_energy_dist = .false. @@ -1197,8 +1238,8 @@ contains ! handle special case of total cross section if (MT == 1) then - xs = xs + mat % density * (ONE-f) * nuc%sigma_t(IE) + & - & f * (nuc%sigma_t(IE+1)) + xs = xs + mat % density * (ONE-f) * nuc%total(IE) + & + & f * (nuc%total(IE+1)) cycle end if @@ -1376,12 +1417,12 @@ contains ! get nuclide grid energy and cross-section value E_i = nuc % energy(IE) - sigma_i = nuc % sigma_t(IE) + sigma_i = nuc % total(IE) ! interpolate as necessary if (E_i < e_grid(k)) then E_i1 = nuc % energy(IE + 1) - sigma_i1 = nuc % sigma_t(IE + 1) + sigma_i1 = nuc % total(IE + 1) r = (e_grid(k) - E_i)/(E_i1 - E_i) val = (ONE - r)*sigma_i + r*sigma_i1 diff --git a/src/cross_section_header.f90 b/src/cross_section_header.f90 index 8b3a2a515c..9368dfabb9 100644 --- a/src/cross_section_header.f90 +++ b/src/cross_section_header.f90 @@ -56,14 +56,20 @@ module cross_section_header type Nuclide character(20) :: name - real(8) :: awr - real(8) :: temp + real(8) :: awr + real(8) :: temp + logical :: fissionable + + ! Energy grid information integer :: n_grid integer, allocatable :: grid_index(:) real(8), allocatable :: energy(:) - real(8), allocatable :: sigma_t(:) - real(8), allocatable :: sigma_a(:) - real(8), allocatable :: sigma_el(:) + + ! Cross sections + real(8), allocatable :: total(:) + real(8), allocatable :: elastic(:) + real(8), allocatable :: fission(:) + real(8), allocatable :: absorption(:) real(8), allocatable :: heating(:) ! Total fission neutron emission @@ -133,4 +139,36 @@ module cross_section_header character(150) :: path end type xsData +!=============================================================================== +! NUCLIDEMICROXS contains cached microscopic cross sections for a +! particular nuclide at the current energy +!=============================================================================== + + type NuclideMicroXS + integer :: index_grid + integer :: index_temp + integer :: last_index_grid + integer :: last_index_temp + real(8) :: interp_factor + real(8) :: total + real(8) :: elastic + real(8) :: absorption + real(8) :: fission + real(8) :: nu_fission + end type NuclideMicroXS + +!=============================================================================== +! MATERIALMACROXS contains cached macroscopic cross sections for the material a +! particle is traveling through +!=============================================================================== + + type MaterialMacroXS + real(8) :: total + real(8) :: scatter + real(8) :: elastic + real(8) :: absorption + real(8) :: fission + real(8) :: nu_fission + end type MaterialMacroXS + end module cross_section_header diff --git a/src/endf.f90 b/src/endf.f90 index 8232e8f5a4..5d172b138b 100644 --- a/src/endf.f90 +++ b/src/endf.f90 @@ -27,7 +27,7 @@ contains string = '(n,2n)' case (N_3N) string = '(n,3n)' - case (FISSION) + case (N_FISSION) string = '(n,fission)' case (N_F) string = '(n,f)' diff --git a/src/fission.f90 b/src/fission.f90 new file mode 100644 index 0000000000..3b577f5e48 --- /dev/null +++ b/src/fission.f90 @@ -0,0 +1,208 @@ +module fission + + use constants + use cross_section_header, only: Nuclide + use error, only: fatal_error + use search, only: binary_search + +contains + +!=============================================================================== +! NU_TOTAL calculates the total number of neutrons emitted per fission for a +! given nuclide and incoming neutron energy +!=============================================================================== + + function nu_total(nuc, E) result(nu) + + type(Nuclide), pointer :: nuc ! nuclide from which to find nu + real(8), intent(in) :: E ! energy of incoming neutron + real(8) :: nu ! number of total neutrons emitted per fission + + integer :: i ! loop index + integer :: j ! index on nu energy grid / precursor group + integer :: loc ! index before start of energies/nu values + integer :: NC ! number of polynomial coefficients + integer :: NR ! number of interpolation regions + integer :: NE ! number of energies tabulated + real(8) :: c ! polynomial coefficient + real(8) :: f ! interpolation factor + character(MAX_LINE_LEN) :: msg ! error message + + if (nuc % nu_t_type == NU_NONE) then + msg = "No neutron emission data for table: " // nuc % name + call fatal_error(msg) + elseif (nuc % nu_t_type == NU_POLYNOMIAL) then + ! determine number of coefficients + NC = int(nuc % nu_t_data(1)) + + ! sum up polynomial in energy + nu = ZERO + do i = 0, NC - 1 + c = nuc % nu_t_data(i+2) + nu = nu + c * E**i + end do + elseif (nuc % nu_t_type == NU_TABULAR) then + ! determine number of interpolation regions -- as far as I can tell, no + ! nu data has multiple interpolation regions. Furthermore, it seems all + ! are lin-lin. + NR = int(nuc % nu_t_data(1)) + if (NR /= 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to determine total nu." + call fatal_error(msg) + end if + + ! determine number of energies + loc = 2 + 2*NR + NE = int(nuc % nu_t_data(loc)) + + ! do binary search over tabuled energies to determine appropriate index + ! and interpolation factor + if (E < nuc % nu_t_data(loc+1)) then + j = 1 + f = ZERO + elseif (E > nuc % nu_t_data(loc+NE)) then + j = NE - 1 + f = ONE + else + j = binary_search(nuc % nu_t_data(loc+1), NE, E) + f = (E - nuc % nu_t_data(loc+j)) / & + & (nuc % nu_t_data(loc+j+1) - nuc % nu_t_data(loc+j)) + end if + + ! determine nu total + loc = loc + NE + nu = nuc % nu_t_data(loc+j) + f * & + & (nuc % nu_t_data(loc+j+1) - nuc % nu_t_data(loc+j)) + end if + + end function nu_total + +!=============================================================================== +! NU_PROMPT calculates the total number of prompt neutrons emitted per fission +! for a given nuclide and incoming neutron energy +!=============================================================================== + + function nu_prompt(nuc, E) result(nu) + + type(Nuclide), pointer :: nuc ! nuclide from which to find nu + real(8), intent(in) :: E ! energy of incoming neutron + real(8) :: nu ! number of prompt neutrons emitted per fission + + integer :: i ! loop index + integer :: j ! index on nu energy grid / precursor group + integer :: loc ! index before start of energies/nu values + integer :: NC ! number of polynomial coefficients + integer :: NR ! number of interpolation regions + integer :: NE ! number of energies tabulated + real(8) :: c ! polynomial coefficient + real(8) :: f ! interpolation factor + character(MAX_LINE_LEN) :: msg ! error message + + if (nuc % nu_p_type == NU_NONE) then + ! since no prompt or delayed data is present, this means all neutron + ! emission is prompt -- WARNING: This currently returns zero. The calling + ! routine needs to know this situation is occurring since we don't want + ! to call nu_total unnecessarily if it's already been called + nu = ZERO + elseif (nuc % nu_p_type == NU_POLYNOMIAL) then + ! determine number of coefficients + NC = int(nuc % nu_p_data(1)) + + ! sum up polynomial in energy + nu = ZERO + do i = 0, NC - 1 + c = nuc % nu_p_data(i+2) + nu = nu + c * E**i + end do + elseif (nuc % nu_p_type == NU_TABULAR) then + ! determine number of interpolation regions + NR = int(nuc % nu_p_data(1)) + if (NR /= 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to determine prompt nu." + call fatal_error(msg) + end if + + ! determine number of energies + loc = 2 + 2*NR + NE = int(nuc % nu_p_data(loc)) + + ! do binary search over tabuled energies to determine appropriate index + ! and interpolation factor + if (E < nuc % nu_p_data(loc+1)) then + j = 1 + f = ZERO + elseif (E > nuc % nu_p_data(loc+NE)) then + j = NE - 1 + f = ONE + else + j = binary_search(nuc % nu_p_data(loc+1), NE, E) + f = (E - nuc % nu_p_data(loc+j)) / & + & (nuc % nu_p_data(loc+j+1) - nuc % nu_p_data(loc+j)) + end if + + ! determine nu total + loc = loc + NE + nu = nuc % nu_p_data(loc+j) + f * & + & (nuc % nu_p_data(loc+j+1) - nuc % nu_p_data(loc+j)) + end if + + end function nu_prompt + +!=============================================================================== +! NU_DELAYED calculates the total number of delayed neutrons emitted per fission +! for a given nuclide and incoming neutron energy +!=============================================================================== + + function nu_delayed(nuc, E) result(nu) + + type(Nuclide), pointer :: nuc ! nuclide from which to find nu + real(8), intent(in) :: E ! energy of incoming neutron + real(8) :: nu ! number of delayed neutrons emitted per fission + + integer :: j ! index on nu energy grid / precursor group + integer :: loc ! index before start of energies/nu values + integer :: NR ! number of interpolation regions + integer :: NE ! number of energies tabulated + real(8) :: f ! interpolation factor + character(MAX_LINE_LEN) :: msg ! error message + + if (nuc % nu_d_type == NU_NONE) then + nu = ZERO + elseif (nuc % nu_d_type == NU_TABULAR) then + ! determine number of interpolation regions + NR = int(nuc % nu_d_data(1)) + if (NR /= 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to determine delayed nu." + call fatal_error(msg) + end if + + ! determine number of energies + loc = 2 + 2*NR + NE = int(nuc % nu_d_data(loc)) + + ! do binary search over tabuled energies to determine appropriate index + ! and interpolation factor + if (E < nuc % nu_d_data(loc+1)) then + j = 1 + f = ZERO + elseif (E > nuc % nu_d_data(loc+NE)) then + j = NE - 1 + f = ONE + else + j = binary_search(nuc % nu_d_data(loc+1), NE, E) + f = (E - nuc % nu_d_data(loc+j)) / & + & (nuc % nu_d_data(loc+j+1) - nuc % nu_d_data(loc+j)) + end if + + ! determine nu total + loc = loc + NE + nu = nuc % nu_d_data(loc+j) + f * & + & (nuc % nu_d_data(loc+j+1) - nuc % nu_d_data(loc+j)) + end if + + end function nu_delayed + +end module fission diff --git a/src/geometry.f90 b/src/geometry.f90 index 6a5d40a71e..cfd9585a56 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -121,8 +121,9 @@ contains cUniverse => univ ! set particle attributes - p%cell = univ % cells(i) - p%universe = dict_get_key(universe_dict, univ % uid) + p % cell = univ % cells(i) + p % universe = dict_get_key(universe_dict, univ % uid) + p % material = c % material exit elseif (c % type == CELL_FILL) then lower_univ => universes(c % fill) @@ -337,7 +338,8 @@ contains end if else ! set current pointers - p%cell = index_cell + p % cell = index_cell + p % material = c % material cCell => c cMaterial => materials(cCell%material) end if @@ -360,7 +362,8 @@ contains end if else ! set current pointers - p%cell = index_cell + p % cell = index_cell + p % material = c % material cCell => c cMaterial => materials(cCell%material) end if @@ -375,7 +378,8 @@ contains do i = 1, size(cells) c => cells(i) if (cell_contains(c, p)) then - p%cell = i + p % cell = i + p % material = c % material cCell => c cMaterial => materials(cCell%material) return diff --git a/src/global.f90 b/src/global.f90 index 2535363cd3..adcbb703cd 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -2,7 +2,8 @@ module global use bank_header, only: Bank use constants - use cross_section_header, only: Nuclide, SAB_Table, xsData + use cross_section_header, only: Nuclide, SAB_Table, xsData, NuclideMicroXS, & + MaterialMacroXS use datatypes_header, only: DictionaryII, DictionaryCI use geometry_header, only: Cell, Universe, Lattice, Surface use material_header, only: Material @@ -51,6 +52,10 @@ module global integer :: n_nuclides_total integer :: n_sab_tables + ! Cross section caches + type(NuclideMicroXS), allocatable :: micro_xs(:) + type(MaterialMacroXS) :: material_xs + ! Current cell, surface, material type(Cell), pointer :: cCell type(Universe), pointer :: cUniverse diff --git a/src/particle_header.f90 b/src/particle_header.f90 index fc2650eb89..e1eb3fb26f 100644 --- a/src/particle_header.f90 +++ b/src/particle_header.f90 @@ -8,23 +8,25 @@ module particle_header !=============================================================================== type Particle - integer(8) :: uid ! Unique ID - integer :: type ! Particle type (n, p, e, etc) - real(8) :: xyz(3) ! location - real(8) :: xyz_local(3) ! local location (after transformations) - real(8) :: uvw(3) ! directional cosines - real(8) :: wgt ! particle weight - real(8) :: E ! energy - integer :: IE ! index on energy grid - real(8) :: interp ! interpolation factor for energy grid - integer :: cell ! current cell - integer :: universe ! current universe - integer :: lattice ! current lattice - integer :: surface ! current surface - integer :: index_x ! lattice index for x direction - integer :: index_y ! lattice index for y direction - logical :: alive ! is particle alive? - integer :: n_coll ! # of collisions + integer(8) :: uid ! Unique ID + integer :: type ! Particle type (n, p, e, etc) + real(8) :: xyz(3) ! location + real(8) :: xyz_local(3) ! local location (after transformations) + real(8) :: uvw(3) ! directional cosines + real(8) :: wgt ! particle weight + real(8) :: E ! energy + integer :: IE ! index on energy grid + real(8) :: interp ! interpolation factor for energy grid + integer :: cell ! current cell + integer :: universe ! current universe + integer :: lattice ! current lattice + integer :: surface ! current surface + integer :: material ! current material + integer :: last_material ! last material + integer :: index_x ! lattice index for x direction + integer :: index_y ! lattice index for y direction + logical :: alive ! is particle alive? + integer :: n_coll ! # of collisions end type Particle end module particle_header diff --git a/src/physics.f90 b/src/physics.f90 index d3b5990434..b376be287d 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -4,6 +4,7 @@ module physics use cross_section_header, only: Nuclide, Reaction, DistEnergy use endf, only: reaction_name use error, only: fatal_error, warning + use fission, only: nu_total, nu_prompt, nu_delayed use geometry, only: find_cell, dist_to_boundary, cross_surface, & cross_lattice use geometry_header, only: Universe, BASE_UNIVERSE @@ -29,12 +30,9 @@ contains integer :: surf ! surface which particle is on integer :: last_cell ! most recent cell particle was in - integer :: IE ! index on energy grid real(8) :: d_to_boundary ! distance to nearest boundary real(8) :: d_to_collision ! sampled distance to collision real(8) :: distance ! distance particle travels - real(8) :: Sigma ! total cross-section - real(8) :: f ! interpolation factor logical :: found_cell ! found cell which particle is in? logical :: in_lattice ! is surface crossing in lattice? character(MAX_LINE_LEN) :: msg ! output/error message @@ -46,8 +44,8 @@ contains ! if particle couldn't be located, print error if (.not. found_cell) then - write(msg, 100) "Could not locate cell for particle at: ", p % xyz -100 format (A,3ES11.3) + write(msg, '(A,3ES11.3)') & + "Could not locate cell for particle at: ", p % xyz call fatal_error(msg) end if end if @@ -63,22 +61,16 @@ contains end if ! find energy index, interpolation factor - call find_energy_index(p) - do while (p % alive) - ! Copy energy index and interpolation factor - IE = p % IE - f = p % interp - - ! Determine material total cross-section - Sigma = (1-f)*cMaterial%total_xs(IE) + f*cMaterial%total_xs(IE+1) + ! Calculate microscopic and macroscopic cross sections + call calculate_xs(p) ! Find the distance to the nearest boundary call dist_to_boundary(p, d_to_boundary, surf, in_lattice) ! Sample a distance to collision - d_to_collision = -log(rang()) / Sigma + d_to_collision = -log(rang()) / material_xs % total ! Select smaller of the two distances distance = min(d_to_boundary, d_to_collision) @@ -107,6 +99,131 @@ contains end subroutine transport +!=============================================================================== +! CALCULATE_XS determines the macroscopic cross sections for the material the +! particle is currently traveling through. +!=============================================================================== + + subroutine calculate_xs(p) + + type(Particle), pointer :: p + + integer :: i ! loop index over nuclides + integer :: index_nuclide ! index into nuclides array + real(8) :: atom_density ! atom density of a nuclide + type(Material), pointer :: mat ! current material + + ! If the material is the same as the last material and the energy of the + ! particle hasn't changed, we don't need to lookup cross sections again. + + if (p % material == p % last_material) return + + ! Set all material macroscopic cross sections to zero + material_xs % total = ZERO + material_xs % scatter = ZERO + material_xs % absorption = ZERO + material_xs % fission = ZERO + material_xs % nu_fission = ZERO + + mat => materials(p % material) + + ! Find energy index on unionized grid + call find_energy_index(p) + + ! Add contribution from each nuclide in material + do i = 1, mat % n_nuclides + ! Determine microscopic cross sections for this nuclide + index_nuclide = mat % nuclide(i) + call calculate_nuclide_xs(p, index_nuclide) + + ! Copy atom density of nuclide in material + atom_density = mat % atom_density(i) + + ! Add contributions to material macroscopic total cross section + material_xs % total = material_xs % total + & + atom_density * micro_xs(index_nuclide) % total + + ! Add contributions to material macroscopic scattering cross section + material_xs % elastic = material_xs % elastic + & + atom_density * micro_xs(index_nuclide) % elastic + + ! Add contributions to material macroscopic absorption cross section + material_xs % absorption = material_xs % absorption + & + atom_density * micro_xs(index_nuclide) % absorption + + ! Add contributions to material macroscopic fission cross section + material_xs % fission = material_xs % fission + & + atom_density * micro_xs(index_nuclide) % fission + + ! Add contributions to material macroscopic nu-fission cross section + material_xs % nu_fission = material_xs % nu_fission + & + atom_density * micro_xs(index_nuclide) % nu_fission + end do + + end subroutine calculate_xs + +!=============================================================================== +! CALCULATE_NUCLIDE_XS determines microscopic cross sections for a nuclide of a +! given index in the nuclides array at the energy of the given particle +!=============================================================================== + + subroutine calculate_nuclide_xs(p, index_nuclide) + + type(Particle), pointer :: p + integer, intent(in) :: index_nuclide ! index into nuclides array + + integer :: i ! index into nuclides array + integer :: IE ! index on nuclide energy grid + real(8) :: f ! interpolation factor on nuclide energy grid + real(8) :: nu_t ! total number of neutrons emitted per fission + type(Nuclide), pointer :: nuc ! pointer to nuclide cross section table + + ! Copy index of nuclide + i = index_nuclide + + ! Set pointer to nuclide + nuc => nuclides(i) + + ! TODO: Check if last energy/temp combination is same as current. If so, we + ! can return. + + ! TODO: If not using unionized energy grid, we need to find the index on the + ! nuclide energy grid using lethargy mapping or whatever other technique + + ! search nuclide energy grid + IE = nuc % grid_index(p % IE) + f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE)) + + micro_xs(i) % index_grid = IE + micro_xs(i) % interp_factor = f + + ! Initialize nuclide cross-sections to zero + micro_xs(i) % nu_fission = ZERO + + ! Calculate microscopic nuclide total cross section + micro_xs(i) % total = & + (ONE-f) * nuc % total(IE) + f * nuc % total(IE+1) + + ! Calculate microscopic nuclide total cross section + micro_xs(i) % elastic = & + (ONE-f) * nuc % elastic(IE) + f * nuc % elastic(IE+1) + + ! Calculate microscopic nuclide absorption cross section + micro_xs(i) % absorption = & + (ONE-f) * nuc % absorption(IE) + f * nuc % absorption(IE+1) + + if (nuc % fissionable) then + ! Calculate microscopic nuclide total cross section + micro_xs(i) % fission = & + (ONE-f) * nuc % fission(IE) + f * nuc % fission(IE+1) + + ! Calculate microscopic nuclide nu-fission cross section + nu_t = nu_total(nuc, p % E) + micro_xs(i) % nu_fission = nu_t * micro_xs(i) % fission + end if + + end subroutine calculate_nuclide_xs + !=============================================================================== ! FIND_ENERGY_INDEX determines the index on the union energy grid and the ! interpolation factor for a particle at a certain energy @@ -153,63 +270,47 @@ contains type(Particle), pointer :: p integer :: i ! index over nuclides in a material - integer :: n_nuclides ! number of nuclides in material + integer :: index_nuclide integer :: IE ! index on nuclide energy grid real(8) :: f ! interpolation factor real(8) :: sigma ! microscopic total xs for nuclide real(8) :: total ! total macroscopic xs for material - real(8) :: density_i ! atom density of nuclide real(8) :: prob ! cumulative probability - real(8) :: r1 ! random number + real(8) :: cutoff ! random number real(8) :: flux ! collision estimator of flux - real(8), allocatable :: Sigma_rxn(:) ! macroscopic xs for each nuclide - character(MAX_LINE_LEN) :: msg ! output/error message + real(8) :: atom_density ! atom density of nuclide in atom/b-cm + character(MAX_LINE_LEN) :: msg ! output/error message type(Material), pointer :: mat type(Nuclide), pointer :: nuc type(Reaction), pointer :: rxn - mat => cMaterial - - ! calculate total cross-section for each nuclide at current energy in order - ! to create discrete pdf for sampling nuclide - n_nuclides = mat % n_nuclides - allocate(Sigma_rxn(n_nuclides)) - do i = 1, n_nuclides - nuc => nuclides(mat % nuclide(i)) - - ! determine nuclide atom density - density_i = mat % atom_density(i) - - ! search nuclide energy grid - IE = nuc%grid_index(p % IE) - f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE)) - - ! calculate nuclide macroscopic cross-section - Sigma_rxn(i) = density_i*((ONE-f)*nuc%Sigma_t(IE) + & - & f*(nuc%Sigma_t(IE+1))) - end do - total = sum(Sigma_rxn) + mat => materials(p % material) ! sample nuclide - r1 = rang() + i = 0 + total = material_xs % total + cutoff = rang() * total prob = ZERO - do i = 1, n_nuclides - prob = prob + Sigma_rxn(i) / total - if (r1 < prob) exit + do while (prob < cutoff) + i = i + 1 + if (i > mat % n_nuclides) then + msg = "Did not sample any nuclide during collision." + call fatal_error(msg) + end if + + index_nuclide = mat % nuclide(i) + atom_density = mat % atom_density(i) + sigma = atom_density * micro_xs(index_nuclide) % total + prob = prob + sigma end do ! Get table, total xs, interpolation factor - density_i = mat % atom_density(i) - nuc => nuclides(mat%nuclide(i)) - sigma = Sigma_rxn(i) / density_i + nuc => nuclides(index_nuclide) IE = nuc%grid_index(p % IE) f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE)) - ! free memory - deallocate(Sigma_rxn) - ! sample reaction channel - r1 = rang()*sigma + cutoff = rang() * sigma / atom_density prob = ZERO do i = 1, nuc % n_reaction rxn => nuc % reactions(i) @@ -225,7 +326,7 @@ contains ! add to cumulative probability prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) & & + f*(rxn%sigma(IE-rxn%IE+2))) - if (r1 < prob) exit + if (cutoff < prob) exit end do if (verbosity >= 10) then @@ -243,8 +344,8 @@ contains case (N_NA, N_N3A, N_NP, N_N2A, N_ND, N_NT, N_N3HE, & & N_NT2A, N_N2P, N_NPA, N_N1 : N_NC) call inelastic_scatter(p, nuc, rxn) - case (FISSION, N_F, N_NF, N_2NF, N_3NF) - call n_fission(p, nuc, rxn) + case (N_FISSION, N_F, N_NF, N_2NF, N_3NF) + call n_fission_event(p, nuc, rxn) case (N_GAMMA : N_DA) call n_absorption(p) case default @@ -350,7 +451,7 @@ contains ! with implicit absorption, namely sampling of the number of neutrons! !=============================================================================== - subroutine n_fission(p, nuc, rxn) + subroutine n_fission_event(p, nuc, rxn) type(Particle), pointer :: p type(Nuclide), pointer :: nuc @@ -360,18 +461,16 @@ contains integer :: j ! index on nu energy grid / precursor group integer :: k ! index on precursor yield grid integer :: loc ! index before start of energies/nu values - integer :: NC ! number of polynomial coefficients integer :: NR ! number of interpolation regions integer :: NE ! number of energies tabulated integer :: nu ! actual number of neutrons produced integer :: law ! energy distribution law real(8) :: E ! incoming energy of neutron real(8) :: E_out ! outgoing energy of fission neutron - real(8) :: c ! polynomial coefficient real(8) :: f ! interpolation factor - real(8) :: nu_total ! total nu - real(8) :: nu_prompt ! prompt nu - real(8) :: nu_delay ! delayed nu + real(8) :: nu_t ! total nu + real(8) :: nu_p ! prompt nu + real(8) :: nu_d ! delayed nu real(8) :: mu ! fission neutron angular cosine real(8) :: phi ! fission neutron azimuthal angle real(8) :: beta ! delayed neutron fraction @@ -383,152 +482,30 @@ contains ! copy energy of neutron E = p % E - ! ========================================================================== - ! DETERMINE TOTAL NU - if (nuc % nu_t_type == NU_NONE) then - msg = "No neutron emission data for table: " // nuc % name - call fatal_error(msg) - elseif (nuc % nu_t_type == NU_POLYNOMIAL) then - ! determine number of coefficients - NC = int(nuc % nu_t_data(1)) + ! Determine total nu + nu_t = nu_total(nuc, E) - ! sum up polynomial in energy - nu_total = ZERO - do i = 0, NC - 1 - c = nuc % nu_t_data(i+2) - nu_total = nu_total + c * E**i - end do - elseif (nuc % nu_t_type == NU_TABULAR) then - ! determine number of interpolation regions -- as far as I can tell, no - ! nu data has multiple interpolation regions. Furthermore, it seems all - ! are lin-lin. - NR = int(nuc % nu_t_data(1)) - if (NR /= 0) then - msg = "Multiple interpolation regions not supported while & - &attempting to determine total nu." - call fatal_error(msg) - end if - - ! determine number of energies - loc = 2 + 2*NR - NE = int(nuc % nu_t_data(loc)) - - ! do binary search over tabuled energies to determine appropriate index - ! and interpolation factor - if (E < nuc % nu_t_data(loc+1)) then - j = 1 - f = ZERO - elseif (E > nuc % nu_t_data(loc+NE)) then - j = NE - 1 - f = ONE - else - j = binary_search(nuc % nu_t_data(loc+1), NE, E) - f = (E - nuc % nu_t_data(loc+j)) / & - & (nuc % nu_t_data(loc+j+1) - nuc % nu_t_data(loc+j)) - end if - - ! determine nu total - loc = loc + NE - nu_total = nuc % nu_t_data(loc+j) + f * & - & (nuc % nu_t_data(loc+j+1) - nuc % nu_t_data(loc+j)) + ! Determine prompt nu + if (nuc % nu_p_type == NU_NONE) then + nu_p = nu_t + else + nu_p = nu_prompt(nuc, E) end if - ! ========================================================================== - ! DETERMINE PROMPT NU - if (nuc % nu_p_type == NU_NONE) then - ! since no prompt or delayed data is present, this means all neutron - ! emission is prompt - nu_prompt = nu_total - elseif (nuc % nu_p_type == NU_POLYNOMIAL) then - ! determine number of coefficients - NC = int(nuc % nu_p_data(1)) + ! Determine delayed nu + nu_d = nu_delayed(nuc, E) - ! sum up polynomial in energy - nu_prompt = ZERO - do i = 0, NC - 1 - c = nuc % nu_p_data(i+2) - nu_prompt = nu_prompt + c * E**i - end do - elseif (nuc % nu_p_type == NU_TABULAR) then - ! determine number of interpolation regions - NR = int(nuc % nu_p_data(1)) - if (NR /= 0) then - msg = "Multiple interpolation regions not supported while & - &attempting to determine prompt nu." - call fatal_error(msg) - end if - - ! determine number of energies - loc = 2 + 2*NR - NE = int(nuc % nu_p_data(loc)) - - ! do binary search over tabuled energies to determine appropriate index - ! and interpolation factor - if (E < nuc % nu_p_data(loc+1)) then - j = 1 - f = ZERO - elseif (E > nuc % nu_p_data(loc+NE)) then - j = NE - 1 - f = ONE - else - j = binary_search(nuc % nu_p_data(loc+1), NE, E) - f = (E - nuc % nu_p_data(loc+j)) / & - & (nuc % nu_p_data(loc+j+1) - nuc % nu_p_data(loc+j)) - end if - - ! determine nu total - loc = loc + NE - nu_prompt = nuc % nu_p_data(loc+j) + f * & - & (nuc % nu_p_data(loc+j+1) - nuc % nu_p_data(loc+j)) - end if - - ! ========================================================================== - ! DETERMINE DELAYED NU - if (nuc % nu_d_type == NU_NONE) then - nu_delay = ZERO - elseif (nuc % nu_d_type == NU_TABULAR) then - ! determine number of interpolation regions - NR = int(nuc % nu_d_data(1)) - if (NR /= 0) then - msg = "Multiple interpolation regions not supported while & - &attempting to determine delayed nu." - call fatal_error(msg) - end if - - ! determine number of energies - loc = 2 + 2*NR - NE = int(nuc % nu_d_data(loc)) - - ! do binary search over tabuled energies to determine appropriate index - ! and interpolation factor - if (E < nuc % nu_d_data(loc+1)) then - j = 1 - f = ZERO - elseif (E > nuc % nu_d_data(loc+NE)) then - j = NE - 1 - f = ONE - else - j = binary_search(nuc % nu_d_data(loc+1), NE, E) - f = (E - nuc % nu_d_data(loc+j)) / & - & (nuc % nu_d_data(loc+j+1) - nuc % nu_d_data(loc+j)) - end if - - ! determine nu total - loc = loc + NE - nu_delay = nuc % nu_d_data(loc+j) + f * & - & (nuc % nu_d_data(loc+j+1) - nuc % nu_d_data(loc+j)) - end if - - beta = nu_delay / nu_total + ! Determine delayed neutron fraction + beta = nu_d / nu_t ! TODO: Heat generation from fission ! Sample number of neutrons produced - nu_total = p%wgt / keff * nu_total - if (rang() > nu_total - int(nu_total)) then - nu = int(nu_total) + nu_t = p%wgt / keff * nu_t + if (rang() > nu_t - int(nu_t)) then + nu = int(nu_t) else - nu = int(nu_total) + 1 + nu = int(nu_t) + 1 end if ! Bank source neutrons @@ -633,7 +610,7 @@ contains ! kill original neutron p % alive = .false. - end subroutine n_fission + end subroutine n_fission_event !=============================================================================== ! INELASTIC_SCATTER handles all reactions with a single secondary neutron (other