diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index bf6f93a3c7..d98804c867 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -174,11 +174,15 @@ should be performed. It has the following attributes/sub-elements: ------------------------- The ```` element determines the treatment of the energy grid during -a simulation. The valid options are "nuclide" and "logarithm". Setting this -element to "nuclide" will cause OpenMC to use a nuclide's energy grid when -determining what points to interpolate between for determining cross sections -(i.e. non-unionized energy grid). Setting this element to "logarithm" causes -OpenMC to use a logarithmic mapping technique described in LA-UR-14-24530_. +a simulation. The valid options are "nuclide", "logarithm", and +"material-union". Setting this element to "nuclide" will cause OpenMC to use a +nuclide's energy grid when determining what points to interpolate between for +determining cross sections (i.e. non-unionized energy grid). Setting this +element to "logarithm" causes OpenMC to use a logarithmic mapping technique +described in LA-UR-14-24530_. Setting this element to "material-union" will +cause OpenMC to create energy grids that are unionized material-by-material and +use these grids when determining the energy-cross section pairs to interpolate +cross section values between. *Default*: logarithm diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b3e1a88779..89a4a60e12 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -37,7 +37,7 @@ endif() set(MPI_ENABLED FALSE) set(HDF5_ENABLED FALSE) -if($ENV{FC} MATCHES "mpi.*") +if($ENV{FC} MATCHES "mpi[^/]*$") message("-- Detected MPI wrapper: $ENV{FC}") add_definitions(-DMPI) set(MPI_ENABLED TRUE) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 1f34a2311b..7753f4734a 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -103,7 +103,7 @@ module ace_header ! Energy grid information integer :: n_grid ! # of nuclide grid points - integer, allocatable :: grid_index(:) ! union grid pointers / log grid mapping + integer, allocatable :: grid_index(:) ! log grid mapping indices real(8), allocatable :: energy(:) ! energy values corresponding to xs ! Microscopic cross sections @@ -372,9 +372,6 @@ module ace_header integer :: i ! Loop counter - if (allocated(this % grid_index)) & - deallocate(this % grid_index) - if (allocated(this % energy)) & deallocate(this % energy, this % total, this % elastic, & & this % fission, this % nu_fission, this % absorption) diff --git a/src/constants.F90 b/src/constants.F90 index bac301e36a..85976c00ac 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -368,8 +368,9 @@ module constants ! Energy grid methods integer, parameter :: & - GRID_NUCLIDE = 1, & ! non-unionized energy grid - GRID_LOGARITHM = 2 ! logarithmic mapping + GRID_NUCLIDE = 1, & ! unique energy grid for each nuclide + GRID_MAT_UNION = 2, & ! material union grids with pointers + GRID_LOGARITHM = 3 ! lethargy mapping ! Running modes integer, parameter :: & diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 1672405e8d..dce4d75f32 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -15,6 +15,9 @@ module cross_section implicit none save + integer :: union_grid_index +!$omp threadprivate(union_grid_index) + contains !=============================================================================== @@ -36,11 +39,11 @@ contains !$omp threadprivate(mat) ! Set all material macroscopic cross sections to zero - material_xs % total = ZERO - material_xs % elastic = ZERO - material_xs % absorption = ZERO - material_xs % fission = ZERO - material_xs % nu_fission = ZERO + material_xs % total = ZERO + material_xs % elastic = ZERO + material_xs % absorption = ZERO + material_xs % fission = ZERO + material_xs % nu_fission = ZERO material_xs % kappa_fission = ZERO ! Exit subroutine if material is void @@ -48,6 +51,10 @@ contains mat => materials(p % material) + ! Find energy index on global or material unionized grid + if (grid_method == GRID_MAT_UNION) & + call find_energy_index(p % E, p % material) + ! Determine if this material has S(a,b) tables check_sab = (mat % n_sab > 0) @@ -88,9 +95,9 @@ contains ! Calculate microscopic cross section for this nuclide if (p % E /= micro_xs(i_nuclide) % last_E) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i) else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i) end if ! ======================================================================== @@ -131,24 +138,32 @@ contains ! given index in the nuclides array at the energy of the given particle !=============================================================================== - subroutine calculate_nuclide_xs(i_nuclide, i_sab, E) + subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat) integer, intent(in) :: i_nuclide ! index into nuclides array integer, intent(in) :: i_sab ! index into sab_tables array - real(8), intent(in) :: E ! energy - - integer :: i_grid ! index on nuclide energy grid - integer :: i_low, i_high ! bounding indices from logarithmic mapping - integer :: u ! index into logarithmic mapping array + integer, intent(in) :: i_mat ! index into materials array + integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material + integer :: i_grid ! index on nuclide energy grid + integer :: i_low ! lower logarithmic mapping index + integer :: i_high ! upper logarithmic mapping index + integer :: u ! index into logarithmic mapping array + real(8), intent(in) :: E ! energy real(8) :: f ! interp factor on nuclide energy grid - type(Nuclide), pointer, save :: nuc => null() -!$omp threadprivate(nuc) + type(Nuclide), pointer, save :: nuc => null() + type(Material), pointer, save :: mat => null() +!$omp threadprivate(nuc, mat) - ! Set pointer to nuclide + ! Set pointer to nuclide and material nuc => nuclides(i_nuclide) + mat => materials(i_mat) ! Determine index on nuclide energy grid select case (grid_method) + case (GRID_MAT_UNION) + + i_grid = mat % nuclide_grid_index(i_nuc_mat, union_grid_index) + case (GRID_LOGARITHM) ! Determine the energy grid index using a logarithmic mapping to reduce ! the energy range over which a binary search needs to be performed @@ -173,7 +188,7 @@ contains ! Perform binary search on the nuclide energy grid in order to determine ! which points to interpolate between - if (E < nuc % energy(1)) then + if (E <= nuc % energy(1)) then i_grid = 1 elseif (E > nuc % energy(nuc % n_grid)) then i_grid = nuc % n_grid - 1 @@ -506,6 +521,32 @@ contains end subroutine calculate_urr_xs +!=============================================================================== +! FIND_ENERGY_INDEX determines the index on the union energy grid at a certain +! energy +!=============================================================================== + + subroutine find_energy_index(E, i_mat) + + real(8), intent(in) :: E ! energy of particle + integer, intent(in) :: i_mat ! material index + type(Material), pointer, save :: mat => null() ! pointer to current material +!$omp threadprivate(mat) + + mat => materials(i_mat) + + ! if the energy is outside of energy grid range, set to first or last + ! index. Otherwise, do a binary search through the union energy grid. + if (E <= mat % e_grid(1)) then + union_grid_index = 1 + elseif (E > mat % e_grid(mat % n_grid)) then + union_grid_index = mat % n_grid - 1 + else + union_grid_index = binary_search(mat % e_grid, mat % n_grid, E) + end if + + end subroutine find_energy_index + !=============================================================================== ! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section ! for a given nuclide at the trial relative energy used in resonance scattering diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 234a856c1e..676875a5b4 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -1,6 +1,9 @@ module energy_grid use global + use list_header, only: ListReal + use output, only: write_message + implicit none @@ -10,6 +13,53 @@ module energy_grid contains +!=============================================================================== +! UNIONIZED_GRID creates a unionized energy grid, for the entire problem or for +! each material, composed of the grids from each nuclide in the entire problem, +! or each material, respectively. Right now, the grid for each nuclide is added +! into a linked list one at a time with an effective insertion sort. Could be +! done with a hash for all energy points and then a quicksort at the end (what +! hash function to use?) +!=============================================================================== + + subroutine unionized_grid() + + integer :: i ! index in nuclides array + integer :: j ! index in materials array + type(ListReal), pointer, save :: list => null() + type(Nuclide), pointer, save :: nuc => null() + type(Material), pointer, save :: mat => null() + !$omp threadprivate(list, nuc, mat) + + call write_message("Creating unionized energy grid...", 5) + + ! add grid points for each nuclide in the material + do j = 1, n_materials + mat => materials(j) + do i = 1, mat % n_nuclides + nuc => nuclides(mat % nuclide(i)) + call add_grid_points(list, nuc % energy) + end do + + ! set size of unionized material energy grid + mat % n_grid = list % size() + + ! create allocated array from linked list + allocate(mat % e_grid(mat % n_grid)) + do i = 1, mat % n_grid + mat % e_grid(i) = list % get_item(i) + end do + + ! delete linked list and dictionary + call list % clear() + deallocate(list) + end do + + ! Set pointers to unionized energy grid for each nuclide + call grid_pointers() + + end subroutine unionized_grid + !=============================================================================== ! LOGARITHMIC_GRID determines a logarithmic mapping for energies to bounding ! indices on a nuclide energy grid @@ -59,4 +109,109 @@ contains end subroutine logarithmic_grid +!=============================================================================== +! ADD_GRID_POINTS adds energy points from the 'energy' array into a linked list +! of points already stored from previous arrays. +!=============================================================================== + + subroutine add_grid_points(list, energy) + + type(ListReal), pointer :: list + real(8), intent(in) :: energy(:) + + integer :: i ! index in energy array + integer :: n ! size of energy array + integer :: current ! current index + real(8) :: E ! actual energy value + + i = 1 + n = size(energy) + + ! If the original list is empty, we need to allocate the first element and + ! store first energy point + if (.not. associated(list)) then + allocate(list) + do i = 1, n + call list % append(energy(i)) + end do + return + end if + + ! Set current index to beginning of the list + current = 1 + + do while (i <= n) + E = energy(i) + + ! If we've reached the end of the grid energy list, add the remaining + ! energy points to the end + if (current > list % size()) then + ! Finish remaining energies + do while (i <= n) + call list % append(energy(i)) + i = i + 1 + end do + exit + end if + + if (E < list % get_item(current)) then + + ! Insert new energy in this position + call list % insert(current, E) + + ! Advance index in linked list and in new energy grid + i = i + 1 + current = current + 1 + + elseif (E == list % get_item(current)) then + ! Found the exact same energy, no need to store duplicates so just + ! skip and move to next index + i = i + 1 + current = current + 1 + else + current = current + 1 + end if + + end do + + end subroutine add_grid_points + +!=============================================================================== +! GRID_POINTERS creates an array of pointers (ints) for each nuclide to link +! each point on the nuclide energy grid to one on a unionized energy grid +!=============================================================================== + + subroutine grid_pointers() + + integer :: i ! loop index for nuclides + integer :: j ! loop index for nuclide energy grid + integer :: k ! loop index for materials + integer :: index_e ! index on union energy grid + real(8) :: union_energy ! energy on union grid + real(8) :: energy ! energy on nuclide grid + type(Nuclide), pointer :: nuc => null() + type(Material), pointer :: mat => null() + + do k = 1, n_materials + mat => materials(k) + allocate(mat % nuclide_grid_index(mat % n_nuclides, mat % n_grid)) + do i = 1, mat % n_nuclides + nuc => nuclides(mat % nuclide(i)) + + index_e = 1 + energy = nuc % energy(index_e) + + do j = 1, mat % n_grid + union_energy = mat % e_grid(j) + if (union_energy >= energy .and. index_e < nuc % n_grid) then + index_e = index_e + 1 + energy = nuc % energy(index_e) + end if + mat % nuclide_grid_index(i,j) = index_e - 1 + end do + end do + end do + + end subroutine grid_pointers + end module energy_grid diff --git a/src/global.F90 b/src/global.F90 index 004c43746e..06be1c992f 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -218,6 +218,7 @@ module global type(Timer) :: time_total ! timer for total run type(Timer) :: time_initialize ! timer for initialization type(Timer) :: time_read_xs ! timer for reading cross sections + type(Timer) :: time_unionize ! timer for material xs-energy grid union type(Timer) :: time_bank ! timer for fission bank synchronization type(Timer) :: time_bank_sample ! timer for fission bank sampling type(Timer) :: time_bank_sendrecv ! timer for fission bank SEND/RECV diff --git a/src/initialize.F90 b/src/initialize.F90 index e39331675c..51dbf5ed52 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -4,7 +4,7 @@ module initialize use bank_header, only: Bank use constants use dict_header, only: DictIntInt, ElemKeyValueII - use energy_grid, only: logarithmic_grid, grid_method + use energy_grid, only: logarithmic_grid, grid_method, unionized_grid use error, only: fatal_error, warning use geometry, only: neighbor_lists use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,& @@ -109,10 +109,17 @@ contains ! Create linked lists for multiple instances of the same nuclide call same_nuclide_list() - ! Construct logarithmic energy grid for cross-sections - if (grid_method == GRID_LOGARITHM) then + ! Construct unionized or log energy grid for cross-sections + select case (grid_method) + case (GRID_NUCLIDE) + continue + case (GRID_MAT_UNION) + call time_unionize % start() + call unionized_grid() + call time_unionize % stop() + case (GRID_LOGARITHM) call logarithmic_grid() - end if + end select ! Allocate and setup tally stride, matching_bins, and tally maps call configure_tallies() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index cdc9176046..1b55b5b9b9 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -213,8 +213,11 @@ contains select case (trim(temp_str)) case ('nuclide') grid_method = GRID_NUCLIDE - case ('union') - call fatal_error("Union energy grid is no longer supported.") + case ('material-union', 'union') + grid_method = GRID_MAT_UNION + if (trim(temp_str) == 'union') & + call warning('Energy grids will be unionized by material. Global& + & energy grid unionization is no longer an allowed option.') case ('logarithm', 'logarithmic', 'log') grid_method = GRID_LOGARITHM case default diff --git a/src/material_header.F90 b/src/material_header.F90 index fca735fd2a..6434fa5950 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -13,6 +13,13 @@ module material_header real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm + ! Energy grid information + integer :: n_grid ! # of union material grid points + real(8), allocatable :: e_grid(:) ! union material grid energies + + ! Unionized energy grid information + integer, allocatable :: nuclide_grid_index(:,:) ! nuclide e_grid pointers + ! S(a,b) data references integer :: n_sab = 0 ! number of S(a,b) tables integer, allocatable :: i_sab_nuclides(:) ! index of corresponding nuclide diff --git a/src/output.F90 b/src/output.F90 index 5599d8c56f..bd86835d10 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1612,22 +1612,26 @@ contains ! write global tallies if (n_realizations > 1) then - write(ou,102) "k-effective (Collision)", global_tallies(K_COLLISION) & - % sum, global_tallies(K_COLLISION) % sum_sq - write(ou,102) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) & - % sum, global_tallies(K_TRACKLENGTH) % sum_sq - write(ou,102) "k-effective (Absorption)", global_tallies(K_ABSORPTION) & - % sum, global_tallies(K_ABSORPTION) % sum_sq - if (n_realizations > 3) write(ou,102) "Combined k-effective", k_combined + if (run_mode == MODE_EIGENVALUE) then + write(ou,102) "k-effective (Collision)", global_tallies(K_COLLISION) & + % sum, global_tallies(K_COLLISION) % sum_sq + write(ou,102) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) & + % sum, global_tallies(K_TRACKLENGTH) % sum_sq + write(ou,102) "k-effective (Absorption)", global_tallies(K_ABSORPTION) & + % sum, global_tallies(K_ABSORPTION) % sum_sq + if (n_realizations > 3) write(ou,102) "Combined k-effective", k_combined + end if write(ou,102) "Leakage Fraction", global_tallies(LEAKAGE) % sum, & global_tallies(LEAKAGE) % sum_sq else if (master) call warning("Could not compute uncertainties -- only one & &active batch simulated!") - write(ou,103) "k-effective (Collision)", global_tallies(K_COLLISION) % sum - write(ou,103) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) % sum - write(ou,103) "k-effective (Absorption)", global_tallies(K_ABSORPTION) % sum + if (run_mode == MODE_EIGENVALUE) then + write(ou,103) "k-effective (Collision)", global_tallies(K_COLLISION) % sum + write(ou,103) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) % sum + write(ou,103) "k-effective (Absorption)", global_tallies(K_ABSORPTION) % sum + end if write(ou,103) "Leakage Fraction", global_tallies(LEAKAGE) % sum end if write(ou,*) diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index a7b92e1738..092bf9983c 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -27,7 +27,7 @@ element settings { (element weight_avg { xsd:double } | attribute weight_avg { xsd:double })? }? & - element energy_grid { ( "nuclide" | "log" | "logarithm" | "logarithmic" ) }? & + element energy_grid { ( "nuclide" | "log" | "logarithm" | "logarithmic" | "material-union" | "union" ) }? & element entropy { (element dimension { list { xsd:int+ } } | diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 02dba3688e..10770f80cc 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -106,6 +106,8 @@ log logarithm logarithmic + material-union + union diff --git a/src/tally.F90 b/src/tally.F90 index c96f1ecd6e..ded646d4da 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -26,6 +26,490 @@ module tally contains +!=============================================================================== +! SCORE_GENERAL adds scores to the tally array for the given filter and nuclide. +! This will work for either analog or tracklength tallies. Note that +! atom_density and flux are not used for analog tallies. +!=============================================================================== + + subroutine score_general(p, t, start_index, filter_index, i_nuclide, & + atom_density, flux) + type(Particle), intent(in) :: p + type(TallyObject), pointer, intent(inout) :: t + integer, intent(in) :: start_index + integer, intent(in) :: i_nuclide + integer, intent(in) :: filter_index ! for % results + real(8), intent(in) :: flux ! flux estimate + real(8), intent(in) :: atom_density ! atom/b-cm + + integer :: i ! loop index for scoring bins + integer :: l ! loop index for nuclides in material + integer :: m ! loop index for reactions + integer :: n ! loop index for legendre order + integer :: num_nm ! Number of N,M orders in harmonic + integer :: q ! loop index for scoring bins + integer :: i_nuc ! index in nuclides array (from material) + integer :: i_energy ! index in nuclide energy grid + integer :: score_bin ! scoring bin, e.g. SCORE_FLUX + integer :: score_index ! scoring bin index + real(8) :: atom_density_ ! atom/b-cm + real(8) :: f ! interpolation factor + real(8) :: score ! analog tally score + real(8) :: macro_total ! material macro total xs + real(8) :: macro_scatt ! material macro scatt xs + real(8) :: uvw(3) ! particle direction + type(Material), pointer, save :: mat => null() + type(Reaction), pointer, save :: rxn => null() +!$omp threadprivate(mat, rxn) + + i = 0 + SCORE_LOOP: do q = 1, t % n_user_score_bins + i = i + 1 + + ! determine what type of score bin + score_bin = t % score_bins(i) + + ! determine scoring bin index + score_index = start_index + i + + !######################################################################### + ! Determine appropirate scoring value. + + select case(score_bin) + + + case (SCORE_FLUX, SCORE_FLUX_YN) + if (t % estimator == ESTIMATOR_ANALOG) then + ! All events score to a flux bin. We actually use a collision + ! estimator since there is no way to count 'events' exactly for + ! the flux + if (survival_biasing) then + ! We need to account for the fact that some weight was already + ! absorbed + score = p % last_wgt + p % absorb_wgt + else + score = p % last_wgt + end if + score = score / material_xs % total + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + ! For flux, we need no cross section + score = flux + end if + + + case (SCORE_TOTAL, SCORE_TOTAL_YN) + if (t % estimator == ESTIMATOR_ANALOG) then + ! All events will score to the total reaction rate. We can just + ! use the weight of the particle entering the collision as the + ! score + if (survival_biasing) then + ! We need to account for the fact that some weight was already + ! absorbed + score = p % last_wgt + p % absorb_wgt + else + score = p % last_wgt + end if + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % total * atom_density * flux + else + score = material_xs % total * flux + end if + end if + + + case (SCORE_SCATTER, SCORE_SCATTER_N) + if (t % estimator == ESTIMATOR_ANALOG) then + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP + ! Since only scattering events make it here, again we can use + ! the weight entering the collision as the estimator for the + ! reaction rate + score = p % last_wgt + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + ! Note SCORE_SCATTER_N not available for tracklength. + if (i_nuclide > 0) then + score = (micro_xs(i_nuclide) % total & + - micro_xs(i_nuclide) % absorption) * atom_density * flux + else + score = (material_xs % total - material_xs % absorption) * flux + end if + end if + + + case (SCORE_SCATTER_PN, SCORE_SCATTER_YN) + ! Only analog estimators are available. + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) then + i = i + t % moment_order(i) + cycle SCORE_LOOP + end if + ! Since only scattering events make it here, again we can use + ! the weight entering the collision as the estimator for the + ! reaction rate + score = p % last_wgt + + + case (SCORE_NU_SCATTER, SCORE_NU_SCATTER_N) + ! Only analog estimators are available. + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP + ! For scattering production, we need to use the post-collision + ! weight as the estimate for the number of neutrons exiting a + ! reaction with neutrons in the exit channel + score = p % wgt + + + case (SCORE_NU_SCATTER_PN, SCORE_NU_SCATTER_YN) + ! Only analog estimators are available. + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) then + i = i + t % moment_order(i) + cycle SCORE_LOOP + end if + ! For scattering production, we need to use the post-collision + ! weight as the estimate for the number of neutrons exiting a + ! reaction with neutrons in the exit channel + score = p % wgt + + + case (SCORE_TRANSPORT) + ! Only analog estimators are available. + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP + ! get material macros + macro_total = material_xs % total + macro_scatt = material_xs % total - material_xs % absorption + ! Score total rate - p1 scatter rate Note estimator needs to be + ! adjusted since tallying is only occuring when a scatter has + ! happened. Effectively this means multiplying the estimator by + ! total/scatter macro + score = (macro_total - p % mu * macro_scatt) * (ONE / macro_scatt) + + + case (SCORE_N_1N) + ! Only analog estimators are available. + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP + ! Skip any events where weight of particle changed + if (p % wgt /= p % last_wgt) cycle SCORE_LOOP + ! All events that reach this point are (n,1n) reactions + score = p % last_wgt + + + case (SCORE_ABSORPTION) + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing) then + ! No absorption events actually occur if survival biasing is on -- + ! just use weight absorbed in survival biasing + score = p % absorb_wgt + else + ! Skip any event where the particle wasn't absorbed + if (p % event == EVENT_SCATTER) cycle SCORE_LOOP + ! All fission and absorption events will contribute here, so we + ! can just use the particle's weight entering the collision + score = p % last_wgt + end if + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % absorption * atom_density * flux + else + score = material_xs % absorption * flux + end if + end if + + + case (SCORE_FISSION) + 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 + if (micro_xs(p % event_nuclide) % absorption > ZERO) then + score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption + else + score = ZERO + end if + 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 reaction rate + score = p % last_wgt * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption + end if + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % fission * atom_density * flux + else + score = material_xs % fission * flux + end if + end if + + + case (SCORE_NU_FISSION) + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing .or. p % fission) then + if (t % find_filter(FILTER_ENERGYOUT) > 0) then + ! Normally, we only need to make contributions to one scoring + ! bin. However, in the case of fission, since multiple fission + ! neutrons were emitted with different energies, multiple + ! outgoing energy bins may have been scored to. The following + ! logic treats this special case and results to multiple bins + call score_fission_eout(p, t, score_index) + cycle SCORE_LOOP + end if + end if + if (survival_biasing) then + ! No fission events occur if survival biasing is on -- need to + ! calculate fraction of absorptions that would have resulted in + ! nu-fission + if (micro_xs(p % event_nuclide) % absorption > ZERO) then + score = p % absorb_wgt * micro_xs(p % event_nuclide) % & + nu_fission / micro_xs(p % event_nuclide) % absorption + else + score = ZERO + end if + else + ! Skip any non-fission events + if (.not. p % fission) cycle SCORE_LOOP + ! If there is no outgoing energy filter, than we only need to + ! score to one bin. For the score to be 'analog', we need to + ! score the number of particles that were banked in the fission + ! bank. Since this was weighted by 1/keff, we multiply by keff + ! to get the proper score. + score = keff * p % wgt_bank + end if + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % nu_fission * atom_density * flux + else + score = material_xs % nu_fission * flux + end if + end if + + + case (SCORE_KAPPA_FISSION) + 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 scale by kappa-fission + if (micro_xs(p % event_nuclide) % absorption > ZERO) then + score = p % absorb_wgt * & + micro_xs(p % event_nuclide) % kappa_fission / & + micro_xs(p % event_nuclide) % absorption + else + score = ZERO + end if + 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 + score = p % last_wgt * & + micro_xs(p % event_nuclide) % kappa_fission / & + micro_xs(p % event_nuclide) % absorption + end if + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % kappa_fission * atom_density * flux + else + score = material_xs % kappa_fission * flux + end if + end if + + + case (SCORE_EVENTS) + ! Simply count number of scoring events + score = ONE + + + case default + if (t % estimator == ESTIMATOR_ANALOG) then + ! Any other score is assumed to be a MT number. Thus, we just need + ! to check if it matches the MT number of the event + if (p % event_MT /= score_bin) cycle SCORE_LOOP + score = p % last_wgt + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + ! Any other cross section has to be calculated on-the-fly. For + ! cross sections that are used often (e.g. n2n, ngamma, etc. for + ! depletion), it might make sense to optimize this section or + ! pre-calculate cross sections + if (score_bin > 1) then + ! Set default score + score = ZERO + + if (i_nuclide > 0) then + ! TODO: The following search for the matching reaction could + ! be replaced by adding a dictionary on each Nuclide instance + ! of the form {MT: i_reaction, ...} + REACTION_LOOP: do m = 1, nuclides(i_nuclide) % n_reaction + ! Get pointer to reaction + rxn => nuclides(i_nuclide) % reactions(m) + ! Check if this is the desired MT + if (score_bin == rxn % MT) then + ! Retrieve index on nuclide energy grid and interpolation + ! factor + i_energy = micro_xs(i_nuclide) % index_grid + f = micro_xs(i_nuclide) % interp_factor + if (i_energy >= rxn % threshold) then + score = ((ONE - f) * rxn % sigma(i_energy - & + rxn%threshold + 1) + f * rxn % sigma(i_energy - & + rxn%threshold + 2)) * atom_density * flux + end if + exit REACTION_LOOP + end if + end do REACTION_LOOP + + else + ! Get pointer to current material + mat => materials(p % material) + do l = 1, mat % n_nuclides + ! Get atom density + atom_density_ = mat % atom_density(l) + ! Get index in nuclides array + i_nuc = mat % nuclide(l) + ! TODO: The following search for the matching reaction could + ! be replaced by adding a dictionary on each Nuclide + ! instance of the form {MT: i_reaction, ...} + do m = 1, nuclides(i_nuc) % n_reaction + ! Get pointer to reaction + rxn => nuclides(i_nuc) % reactions(m) + ! Check if this is the desired MT + if (score_bin == rxn % MT) then + ! Retrieve index on nuclide energy grid and interpolation + ! factor + i_energy = micro_xs(i_nuc) % index_grid + f = micro_xs(i_nuc) % interp_factor + if (i_energy >= rxn % threshold) then + score = score + ((ONE - f) * rxn % sigma(i_energy - & + rxn%threshold + 1) + f * rxn % sigma(i_energy - & + rxn%threshold + 2)) * atom_density_ * flux + end if + exit + end if + end do + end do + end if + + else + call fatal_error("Invalid score type on tally " & + // to_str(t % id) // ".") + end if + end if + + + end select + + !######################################################################### + ! Expand score if necessary and add to tally results. + + select case(score_bin) + + + case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) + ! Find the scattering order for a singly requested moment, and + ! store its moment contribution. + if (t % moment_order(i) == 1) then + score = score * p % mu ! avoid function call overhead + else + score = score * calc_pn(t % moment_order(i), p % mu) + endif +!$omp atomic + t % results(score_index, filter_index) % value = & + t % results(score_index, filter_index) % value + score + + + case(SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN) + score_index = score_index - 1 + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(i) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) & + % value = t & + % results(score_index: score_index + num_nm - 1, filter_index)& + % value & + + score * calc_pn(n, p % mu) * calc_rn(n, p % last_uvw) +!$omp end critical + end do + i = i + (t % moment_order(i) + 1)**2 - 1 + + + case(SCORE_FLUX_YN, SCORE_TOTAL_YN) + score_index = score_index - 1 + num_nm = 1 + if (t % estimator == ESTIMATOR_ANALOG) then + uvw = p % last_uvw + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + uvw = p % coord0 % uvw + end if + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(i) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) & + % value = t & + % results(score_index: score_index + num_nm - 1, filter_index)& + % value & + + score * calc_rn(n, uvw) +!$omp end critical + end do + i = i + (t % moment_order(i) + 1)**2 - 1 + + + case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN) + score_index = score_index - 1 + ! Find the scattering order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(i) + ! determine scoring bin index + score_index = score_index + 1 + + ! get the score and tally it +!$omp atomic + t % results(score_index, filter_index) % value = & + t % results(score_index, filter_index) % value & + + score * calc_pn(n, p % mu) + end do + i = i + t % moment_order(i) + + + case default +!$omp atomic + t % results(score_index, filter_index) % value = & + t % results(score_index, filter_index) % value + score + + + end select + end do SCORE_LOOP + end subroutine score_general + !=============================================================================== ! SCORE_ANALOG_TALLY keeps track of how many events occur in a specified cell, ! energy range, etc. Note that since these are "analog" tallies, they are only @@ -38,22 +522,13 @@ contains integer :: i integer :: i_tally - integer :: j ! loop index for scoring bins integer :: k ! loop index for nuclide bins - integer :: n ! loop index for legendre order - integer :: num_nm ! Number of N,M orders in harmonic - integer :: l ! scoring bin loop index, allowing for changing ! position during the loop integer :: filter_index ! single index for single bin - integer :: score_bin ! scoring bin, e.g. SCORE_FLUX integer :: i_nuclide ! index in nuclides array - integer :: score_index ! scoring bin index - real(8) :: score ! analog tally score real(8) :: last_wgt ! pre-collision particle weight real(8) :: wgt ! post-collision particle weight real(8) :: mu ! cosine of angle of collision - real(8) :: macro_total ! material macro total xs - real(8) :: macro_scatt ! material macro scatt xs logical :: found_bin ! scoring bin found? type(TallyObject), pointer, save :: t => null() !$omp threadprivate(t) @@ -124,427 +599,8 @@ contains end if ! Determine score for each bin - j = 0 - SCORE_LOOP: do l = 1, t % n_user_score_bins - j = j + 1 - ! determine what type of score bin - score_bin = t % score_bins(j) - - ! determine scoring bin index - score_index = (k - 1)*t % n_score_bins + j - - select case (score_bin) - case (SCORE_FLUX) - ! All events score to a flux bin. We actually use a collision - ! estimator since there is no way to count 'events' exactly for - ! the flux - - if (survival_biasing) then - ! We need to account for the fact that some weight was already - ! absorbed - score = last_wgt + p % absorb_wgt - else - score = last_wgt - end if - - score = score / material_xs % total - case (SCORE_FLUX_YN) - ! All events score to a flux bin. We actually use a collision - ! estimator since there is no way to count 'events' exactly for - ! the flux - - score_index = score_index - 1 - - ! get the score - if (survival_biasing) then - ! We need to account for the fact that some weight was already - ! absorbed - score = last_wgt + p % absorb_wgt - else - score = last_wgt - end if - - score = score / material_xs % total - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % last_uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TOTAL) - ! All events will score to the total reaction rate. We can just - ! use the weight of the particle entering the collision as the - ! score - - if (survival_biasing) then - ! We need to account for the fact that some weight was already - ! absorbed - score = last_wgt + p % absorb_wgt - else - score = last_wgt - end if - - case (SCORE_TOTAL_YN) - ! All events will score to the total reaction rate. We can just - ! use the weight of the particle entering the collision as the - ! score - - score_index = score_index - 1 - - ! get the score - if (survival_biasing) then - ! We need to account for the fact that some weight was already - ! absorbed - score = last_wgt + p % absorb_wgt - else - score = last_wgt - end if - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % last_uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_SCATTER) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! Since only scattering events make it here, again we can use - ! the weight entering the collision as the estimator for the - ! reaction rate - - score = last_wgt - - case (SCORE_NU_SCATTER) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! For scattering production, we need to use the post-collision - ! weight as the estimate for the number of neutrons exiting a - ! reaction with neutrons in the exit channel - - score = wgt - - case (SCORE_SCATTER_N) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! Find the scattering order for a singly requested moment, and - ! store its moment contribution. - - if (t % moment_order(j) == 1) then - score = last_wgt * mu ! avoid function call overhead - else - score = last_wgt * calc_pn(t % moment_order(j), mu) - endif - - case (SCORE_SCATTER_PN) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) then - j = j + t % moment_order(j) - cycle SCORE_LOOP - end if - score_index = score_index - 1 - ! Find the scattering order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + 1 - ! get the score and tally it - score = last_wgt * calc_pn(n, mu) - -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - end do - j = j + t % moment_order(j) - cycle SCORE_LOOP - - case (SCORE_SCATTER_YN) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) then - j = j + t % moment_order(j) - cycle SCORE_LOOP - end if - score_index = score_index - 1 - - ! Calculate the number of moments from t % moment_order - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - ! get the score of the scattering moment - score = last_wgt * calc_pn(n, mu) - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % last_uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_NU_SCATTER_N) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! Find the scattering order for a singly requested moment, and - ! store its moment contribution. - - if (t % moment_order(j) == 1) then - score = wgt * mu ! avoid function call overhead - else - score = wgt * calc_pn(t % moment_order(j), mu) - endif - - case (SCORE_NU_SCATTER_PN) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) then - j = j + t % moment_order(j) - cycle SCORE_LOOP - end if - score_index = score_index - 1 - ! Find the scattering order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + 1 - ! get the score and tally it - score = wgt * calc_pn(n, mu) - -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - end do - j = j + t % moment_order(j) - cycle SCORE_LOOP - - case (SCORE_NU_SCATTER_YN) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) then - j = j + t % moment_order(j) - cycle SCORE_LOOP - end if - score_index = score_index - 1 - - ! Calculate the number of moments from t % moment_order - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - ! get the score of the scattering moment - score = wgt * calc_pn(n, mu) - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % last_uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TRANSPORT) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! get material macros - macro_total = material_xs % total - macro_scatt = material_xs % total - material_xs % absorption - - ! Score total rate - p1 scatter rate Note estimator needs to be - ! adjusted since tallying is only occuring when a scatter has - ! happened. Effectively this means multiplying the estimator by - ! total/scatter macro - score = (macro_total - mu*macro_scatt)*(ONE/macro_scatt) - - case (SCORE_N_1N) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! Skip any events where weight of particle changed - if (wgt /= last_wgt) cycle SCORE_LOOP - - ! All events that reach this point are (n,1n) reactions - score = last_wgt - - case (SCORE_ABSORPTION) - if (survival_biasing) then - ! No absorption events actually occur if survival biasing is on -- - ! just use weight absorbed in survival biasing - - score = p % absorb_wgt - - else - ! Skip any event where the particle wasn't absorbed - if (p % event == EVENT_SCATTER) cycle SCORE_LOOP - - ! All fission and absorption events will contribute here, so we - ! can just use the particle's weight entering the collision - - score = last_wgt - end if - - case (SCORE_FISSION) - 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 - - if (micro_xs(p % event_nuclide) % absorption > ZERO) then - score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission / & - micro_xs(p % event_nuclide) % absorption - else - score = ZERO - end if - - 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 reaction rate - - score = last_wgt * micro_xs(p % event_nuclide) % fission / & - micro_xs(p % event_nuclide) % absorption - end if - - case (SCORE_NU_FISSION) - if (survival_biasing) then - ! No fission events occur if survival biasing is on -- need to - ! calculate fraction of absorptions that would have resulted in - ! nu-fission - - if (t % find_filter(FILTER_ENERGYOUT) > 0) then - ! Normally, we only need to make contributions to one scoring - ! bin. However, in the case of fission, since multiple fission - ! neutrons were emitted with different energies, multiple - ! outgoing energy bins may have been scored to. The following - ! logic treats this special case and results to multiple bins - - call score_fission_eout(p, t, score_index) - cycle SCORE_LOOP - - else - - if (micro_xs(p % event_nuclide) % absorption > ZERO) then - score = p % absorb_wgt * micro_xs(p % event_nuclide) % & - nu_fission / micro_xs(p % event_nuclide) % absorption - else - score = ZERO - end if - end if - - else - ! Skip any non-fission events - if (.not. p % fission) cycle SCORE_LOOP - - if (t % find_filter(FILTER_ENERGYOUT) > 0) then - ! Normally, we only need to make contributions to one scoring - ! bin. However, in the case of fission, since multiple fission - ! neutrons were emitted with different energies, multiple - ! outgoing energy bins may have been scored to. The following - ! logic treats this special case and results to multiple bins - - call score_fission_eout(p, t, score_index) - cycle SCORE_LOOP - - else - ! If there is no outgoing energy filter, than we only need to - ! score to one bin. For the score to be 'analog', we need to - ! score the number of particles that were banked in the fission - ! bank. Since this was weighted by 1/keff, we multiply by keff - ! to get the proper score. - - score = keff * p % wgt_bank - - end if - end if - - case (SCORE_KAPPA_FISSION) - 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 scale by kappa-fission - - if (micro_xs(p % event_nuclide) % absorption > ZERO) then - score = p % absorb_wgt * & - micro_xs(p % event_nuclide) % kappa_fission / & - micro_xs(p % event_nuclide) % absorption - else - score = ZERO - end if - - 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 - score = last_wgt * & - micro_xs(p % event_nuclide) % kappa_fission / & - micro_xs(p % event_nuclide) % absorption - end if - case (SCORE_EVENTS) - ! Simply count number of scoring events - score = ONE - - case default - ! Any other score is assumed to be a MT number. Thus, we just need - ! to check if it matches the MT number of the event - if (p % event_MT /= score_bin) cycle SCORE_LOOP - - score = last_wgt - - end select - - ! Add score to tally -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - - end do SCORE_LOOP + call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & + i_nuclide, ZERO, ZERO) end do NUCLIDE_LOOP @@ -640,26 +696,14 @@ contains integer :: i_tally integer :: j ! loop index for scoring bins integer :: k ! loop index for nuclide bins - integer :: l ! loop index for nuclides in material - integer :: m ! loop index for reactions - integer :: n ! loop index for legendre order - integer :: num_nm ! Number of N,M orders in harmonic - integer :: q ! loop index for scoring bins integer :: filter_index ! single index for single bin integer :: i_nuclide ! index in nuclides array (from bins) - integer :: i_nuc ! index in nuclides array (from material) - integer :: i_energy ! index in nuclide energy grid - integer :: score_bin ! scoring type, e.g. SCORE_FLUX - integer :: score_index ! scoring bin index - real(8) :: f ! interpolation factor real(8) :: flux ! tracklength estimate of flux - real(8) :: score ! actual score (e.g., flux*xs) real(8) :: atom_density ! atom density of single nuclide in atom/b-cm logical :: found_bin ! scoring bin found? type(TallyObject), pointer, save :: t => null() type(Material), pointer, save :: mat => null() - type(Reaction), pointer, save :: rxn => null() -!$omp threadprivate(t, mat, rxn) +!$omp threadprivate(t, mat) ! Determine track-length estimate of flux flux = p % wgt * distance @@ -731,291 +775,8 @@ contains end if ! Determine score for each bin - j = 0 - SCORE_LOOP: do q = 1, t % n_user_score_bins - j = j + 1 - ! determine what type of score bin - score_bin = t % score_bins(j) - - ! determine scoring bin index - score_index = (k - 1)*t % n_score_bins + j - - if (i_nuclide > 0) then - ! ================================================================ - ! DETERMINE NUCLIDE CROSS SECTION - - select case(score_bin) - case (SCORE_FLUX) - ! For flux, we need no cross section - score = flux - - case (SCORE_FLUX_YN) - score_index = score_index - 1 - - ! For flux, we need no cross section - score = flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TOTAL) - ! Total cross section is pre-calculated - score = micro_xs(i_nuclide) % total * & - atom_density * flux - - case (SCORE_TOTAL_YN) - score_index = score_index - 1 - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_SCATTER) - ! Scattering cross section is pre-calculated - score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * & - atom_density * flux - - case (SCORE_ABSORPTION) - ! Absorption cross section is pre-calculated - score = micro_xs(i_nuclide) % absorption * & - atom_density * flux - - case (SCORE_FISSION) - ! Fission cross section is pre-calculated - score = micro_xs(i_nuclide) % fission * & - atom_density * flux - - case (SCORE_NU_FISSION) - ! Nu-fission cross section is pre-calculated - score = micro_xs(i_nuclide) % nu_fission * & - atom_density * flux - - case (SCORE_KAPPA_FISSION) - score = micro_xs(i_nuclide) % kappa_fission * & - atom_density * flux - - case (SCORE_EVENTS) - ! For number of events, just score unity - score = ONE - - case default - ! Any other cross section has to be calculated on-the-fly. For - ! cross sections that are used often (e.g. n2n, ngamma, etc. for - ! depletion), it might make sense to optimize this section or - ! pre-calculate cross sections - - if (score_bin > 1) then - ! Set default score - score = ZERO - - ! TODO: The following search for the matching reaction could - ! be replaced by adding a dictionary on each Nuclide instance - ! of the form {MT: i_reaction, ...} - - REACTION_LOOP: do m = 1, nuclides(i_nuclide) % n_reaction - ! Get pointer to reaction - rxn => nuclides(i_nuclide) % reactions(m) - - ! Check if this is the desired MT - if (score_bin == rxn % MT) then - ! Retrieve index on nuclide energy grid and interpolation - ! factor - i_energy = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - - if (i_energy >= rxn % threshold) then - score = ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density * flux - end if - - exit REACTION_LOOP - end if - end do REACTION_LOOP - - else - call fatal_error("Invalid score type on tally " & - &// to_str(t % id) // ".") - end if - end select - - else - ! ================================================================ - ! DETERMINE MATERIAL CROSS SECTION - - select case(score_bin) - case (SCORE_FLUX) - ! For flux, we need no cross section - score = flux - - case (SCORE_FLUX_YN) - score_index = score_index - 1 - - ! For flux, we need no cross section - score = flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TOTAL) - ! Total cross section is pre-calculated - score = material_xs % total * flux - - case (SCORE_TOTAL_YN) - score_index = score_index - 1 - - ! Total cross section is pre-calculated - score = material_xs % total * flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_SCATTER) - ! Scattering cross section is pre-calculated - score = (material_xs % total - material_xs % absorption) * flux - - case (SCORE_ABSORPTION) - ! Absorption cross section is pre-calculated - score = material_xs % absorption * flux - - case (SCORE_FISSION) - ! Fission cross section is pre-calculated - score = material_xs % fission * flux - - case (SCORE_NU_FISSION) - ! Nu-fission cross section is pre-calculated - score = material_xs % nu_fission * flux - - case (SCORE_KAPPA_FISSION) - score = material_xs % kappa_fission * flux - - case (SCORE_EVENTS) - ! For number of events, just score unity - score = ONE - - case default - ! Any other cross section has to be calculated on-the-fly. This - ! is somewhat costly since it requires a loop over each nuclide - ! in a material and each reaction in the nuclide - - if (score_bin > 1) then - ! Set default score - score = ZERO - - ! Get pointer to current material - mat => materials(p % material) - - do l = 1, mat % n_nuclides - ! Get atom density - atom_density = mat % atom_density(l) - - ! Get index in nuclides array - i_nuc = mat % nuclide(l) - - ! TODO: The following search for the matching reaction could - ! be replaced by adding a dictionary on each Nuclide - ! instance of the form {MT: i_reaction, ...} - - do m = 1, nuclides(i_nuc) % n_reaction - ! Get pointer to reaction - rxn => nuclides(i_nuc) % reactions(m) - - ! Check if this is the desired MT - if (score_bin == rxn % MT) then - ! Retrieve index on nuclide energy grid and interpolation - ! factor - i_energy = micro_xs(i_nuc) % index_grid - f = micro_xs(i_nuc) % interp_factor - - if (i_energy >= rxn % threshold) then - score = score + ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density * flux - end if - - exit - end if - end do - - end do - - else - call fatal_error("Invalid score type on tally " & - &// to_str(t % id) // ".") - end if - end select - end if - - ! Add score to tally -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - - end do SCORE_LOOP + call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & + i_nuclide, atom_density, flux) end do NUCLIDE_BIN_LOOP end if @@ -1047,22 +808,11 @@ contains integer, intent(in) :: filter_index integer :: i ! loop index for nuclides in material - integer :: j ! loop index for scoring bin types - integer :: m ! loop index for reactions in nuclide - integer :: n ! loop index for legendre order - integer :: num_nm ! Number of N,M orders in harmonic - integer :: q ! loop index for scoring bins integer :: i_nuclide ! index in nuclides array - integer :: score_bin ! type of score, e.g. SCORE_FLUX - integer :: score_index ! scoring bin index - integer :: i_energy ! index in nuclide energy grid - real(8) :: f ! interpolation factor - real(8) :: score ! actual scoring tally value real(8) :: atom_density ! atom density of single nuclide in atom/b-cm type(TallyObject), pointer, save :: t => null() type(Material), pointer, save :: mat => null() - type(Reaction), pointer, save :: rxn => null() -!$omp threadprivate(t, mat, rxn) +!$omp threadprivate(t, mat) ! Get pointer to tally t => tallies(i_tally) @@ -1081,290 +831,21 @@ contains i_nuclide = mat % nuclide(i) atom_density = mat % atom_density(i) - ! Loop over score types for each bin - j = 0 - SCORE_LOOP: do q = 1, t % n_user_score_bins - j = j + 1 - ! determine what type of score bin - score_bin = t % score_bins(j) - - ! Determine scoring bin index based on what the index of the nuclide - ! is in the nuclides array - score_index = (i_nuclide - 1)*t % n_score_bins + j - - ! Determine macroscopic nuclide cross section - select case(score_bin) - case (SCORE_FLUX) - score = flux - - case (SCORE_FLUX_YN) - score_index = score_index - 1 - - ! For flux, we need no cross section - score = flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TOTAL) - score = micro_xs(i_nuclide) % total * atom_density * flux - - case (SCORE_TOTAL_YN) - score_index = score_index - 1 - - score = micro_xs(i_nuclide) % total * atom_density * flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_SCATTER) - score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * atom_density * flux - - case (SCORE_ABSORPTION) - score = micro_xs(i_nuclide) % absorption * atom_density * flux - - case (SCORE_FISSION) - score = micro_xs(i_nuclide) % fission * atom_density * flux - - case (SCORE_NU_FISSION) - score = micro_xs(i_nuclide) % nu_fission * atom_density * flux - - case (SCORE_KAPPA_FISSION) - score = micro_xs(i_nuclide) % kappa_fission * atom_density * flux - - case (SCORE_EVENTS) - score = ONE - - case default - ! Any other cross section has to be calculated on-the-fly. For cross - ! sections that are used often (e.g. n2n, ngamma, etc. for depletion), - ! it might make sense to optimize this section or pre-calculate cross - ! sections - - if (score_bin > 1) then - ! Set default score - score = ZERO - - ! TODO: The following search for the matching reaction could be - ! replaced by adding a dictionary on each Nuclide instance of the - ! form {MT: i_reaction, ...} - - REACTION_LOOP: do m = 1, nuclides(i_nuclide) % n_reaction - ! Get pointer to reaction - rxn => nuclides(i_nuclide) % reactions(m) - - ! Check if this is the desired MT - if (score_bin == rxn % MT) then - ! Retrieve index on nuclide energy grid and interpolation factor - i_energy = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - - if (i_energy >= rxn % threshold) then - score = ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density * flux - end if - - exit REACTION_LOOP - end if - end do REACTION_LOOP - - else - call fatal_error("Invalid score type on tally " & - &// to_str(t % id) // ".") - end if - end select - - ! Add score to tally -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - - end do SCORE_LOOP + ! Determine score for each bin + call score_general(p, t, (i_nuclide-1)*t % n_score_bins, filter_index, & + i_nuclide, atom_density, flux) end do NUCLIDE_LOOP ! ========================================================================== ! SCORE TOTAL MATERIAL REACTION RATES - ! Loop over score types for each bin - j = 0 - MATERIAL_SCORE_LOOP: do q = 1, t % n_user_score_bins - j = j + 1 - ! determine what type of score bin - score_bin = t % score_bins(j) + i_nuclide = -1 + atom_density = ZERO - ! Determine scoring bin index based on what the index of the nuclide - ! is in the nuclides array - score_index = n_nuclides_total*t % n_score_bins + j - - ! Determine macroscopic material cross section - select case(score_bin) - case (SCORE_FLUX) - score = flux - - case (SCORE_FLUX_YN) - score_index = score_index - 1 - - ! For flux, we need no cross section - score = flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle MATERIAL_SCORE_LOOP - - case (SCORE_TOTAL) - score = material_xs % total * flux - - case (SCORE_TOTAL_YN) - score_index = score_index - 1 - - ! Total cross section is pre-calculated - score = material_xs % total * flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle MATERIAL_SCORE_LOOP - - case (SCORE_SCATTER) - score = (material_xs % total - material_xs % absorption) * flux - - case (SCORE_ABSORPTION) - score = material_xs % absorption * flux - - case (SCORE_FISSION) - score = material_xs % fission * flux - - case (SCORE_NU_FISSION) - score = material_xs % nu_fission * flux - - case (SCORE_KAPPA_FISSION) - score = material_xs % kappa_fission * flux - - case (SCORE_EVENTS) - score = ONE - - case default - ! Any other cross section has to be calculated on-the-fly. This is - ! somewhat costly since it requires a loop over each nuclide in a - ! material and each reaction in the nuclide - - if (score_bin > 1) then - ! Set default score - score = ZERO - - ! Get pointer to current material - mat => materials(p % material) - - do i = 1, mat % n_nuclides - ! Get atom density - atom_density = mat % atom_density(i) - - ! Get index in nuclides array - i_nuclide = mat % nuclide(i) - - ! TODO: The following search for the matching reaction could - ! be replaced by adding a dictionary on each Nuclide - ! instance of the form {MT: i_reaction, ...} - - do m = 1, nuclides(i_nuclide) % n_reaction - ! Get pointer to reaction - rxn => nuclides(i_nuclide) % reactions(m) - - ! Check if this is the desired MT - if (score_bin == rxn % MT) then - ! Retrieve index on nuclide energy grid and interpolation - ! factor - i_energy = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - - if (i_energy >= rxn % threshold) then - score = score + ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density * flux - end if - - exit - end if - end do - - end do - - else - call fatal_error("Invalid score type on tally " & - &// to_str(t % id) // ".") - end if - end select - - ! Add score to tally -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - - end do MATERIAL_SCORE_LOOP + ! Determine score for each bin + call score_general(p, t, n_nuclides_total*t % n_score_bins, filter_index, & + i_nuclide, atom_density, flux) end subroutine score_all_nuclides @@ -1384,21 +865,15 @@ contains integer :: j ! loop index for direction integer :: k ! loop index for mesh cell crossings integer :: b ! loop index for nuclide bins - integer :: n ! loop index for legendre order - integer :: num_nm ! Number of N,M orders in harmonic - integer :: q ! loop index for scoring bins integer :: ijk0(3) ! indices of starting coordinates integer :: ijk1(3) ! indices of ending coordinates integer :: ijk_cross(3) ! indices of mesh cell crossed integer :: n_cross ! number of surface crossings integer :: filter_index ! single index for single bin - integer :: score_bin ! scoring bin, e.g. SCORE_FLUX integer :: i_nuclide ! index in nuclides array - integer :: score_index ! scoring bin index integer :: i_filter_mesh ! index of mesh filter in filters array real(8) :: atom_density ! density of individual nuclide in atom/b-cm real(8) :: flux ! tracklength estimate of flux - real(8) :: score ! actual score (e.g., flux*xs) real(8) :: uvw(3) ! cosine of angle of particle real(8) :: xyz0(3) ! starting/intermediate coordinates real(8) :: xyz1(3) ! ending coordinates of particle @@ -1612,179 +1087,8 @@ contains end if ! Determine score for each bin - j = 0 - SCORE_LOOP: do q = 1, t % n_user_score_bins - j = j + 1 - ! determine what type of score bin - score_bin = t % score_bins(j) - - ! Determine scoring bin index - score_index = (b - 1)*t % n_score_bins + j - - if (i_nuclide > 0) then - ! Determine macroscopic nuclide cross section - select case(score_bin) - case (SCORE_FLUX) - score = flux - - case (SCORE_FLUX_YN) - score_index = score_index - 1 - - score = flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TOTAL) - score = micro_xs(i_nuclide) % total * & - atom_density * flux - - case (SCORE_TOTAL_YN) - score_index = score_index - 1 - - ! Total cross section is pre-calculated - score = micro_xs(i_nuclide) % total * & - atom_density * flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_SCATTER) - score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * & - atom_density * flux - case (SCORE_ABSORPTION) - score = micro_xs(i_nuclide) % absorption * & - atom_density * flux - case (SCORE_FISSION) - score = micro_xs(i_nuclide) % fission * & - atom_density * flux - case (SCORE_NU_FISSION) - score = micro_xs(i_nuclide) % nu_fission * & - atom_density * flux - case (SCORE_KAPPA_FISSION) - score = micro_xs(i_nuclide) % kappa_fission * atom_density * flux - case (SCORE_EVENTS) - score = ONE - case default - call fatal_error("Invalid score type on tally " & - &// to_str(t % id) // ".") - end select - - else - ! Determine macroscopic material cross section - select case(score_bin) - case (SCORE_FLUX) - score = flux - - case (SCORE_FLUX_YN) - score_index = score_index - 1 - - score = flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TOTAL) - score = material_xs % total * flux - - case (SCORE_TOTAL_YN) - score_index = score_index - 1 - - ! Total cross section is pre-calculated - score = material_xs % total * flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_SCATTER) - score = (material_xs % total - material_xs % absorption) * flux - case (SCORE_ABSORPTION) - score = material_xs % absorption * flux - case (SCORE_FISSION) - score = material_xs % fission * flux - case (SCORE_NU_FISSION) - score = material_xs % nu_fission * flux - case (SCORE_KAPPA_FISSION) - score = material_xs % kappa_fission * flux - case (SCORE_EVENTS) - score = ONE - case default - call fatal_error("Invalid score type on tally " & - &// to_str(t % id) // ".") - end select - end if - - ! Add score to tally -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - - end do SCORE_LOOP + call score_general(p, t, (b-1)*t % n_score_bins, filter_index, & + i_nuclide, atom_density, flux) end do NUCLIDE_BIN_LOOP end if diff --git a/tests/run_tests.py b/tests/run_tests.py index 531225d9ad..42fb846d88 100755 --- a/tests/run_tests.py +++ b/tests/run_tests.py @@ -179,6 +179,10 @@ class Test(object): # Runs cmake when in non-script mode def run_cmake(self): os.environ['FC'] = self.fc + if self.petsc: + os.environ['PETSC_DIR'] = PETSC_DIR + if self.mpi: + os.environ['MPI_DIR'] = MPI_DIR build_opts = self.build_opts.split() self.cmake += build_opts rc = call(self.cmake) diff --git a/tests/test_score_total_yn/results_true.dat b/tests/test_score_total_yn/results_true.dat index b78df10085..36d7040bd9 100644 --- a/tests/test_score_total_yn/results_true.dat +++ b/tests/test_score_total_yn/results_true.dat @@ -60,6 +60,106 @@ tally 2: 0.000000E+00 0.000000E+00 0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +8.795501E-01 +1.600341E-01 +-1.530661E-02 +5.215322E-04 +1.571095E-02 +9.166057E-04 +7.226989E-03 +3.683499E-04 +-2.456165E-02 +1.571441E-04 +1.387449E-02 +5.600648E-04 +-2.186870E-02 +2.081164E-04 +-1.106814E-02 +1.065144E-04 +-4.579593E-03 +2.008225E-04 +7.573253E-03 +2.493075E-04 +-1.505016E-02 +1.451400E-04 +1.503792E-02 +9.539029E-05 +-1.183668E-02 +1.741393E-04 +4.060912E-03 +5.890591E-05 +1.789747E-02 +8.239749E-05 +-2.168438E-02 +1.555757E-04 +-1.045826E-02 +8.108402E-05 +1.227413E-02 +2.502871E-04 +-2.806122E-02 +1.886120E-04 +-4.918718E-03 +2.129038E-04 +-1.014729E-02 +6.914145E-05 +-5.873760E-03 +2.229738E-04 +6.666510E-03 +1.394147E-04 +-3.245528E-03 +1.550876E-04 +-1.037659E-02 +2.163837E-04 1.517577E+01 4.747271E+01 -2.463504E-01 @@ -110,6 +210,56 @@ tally 2: 2.383142E-02 -1.463402E-01 1.526988E-02 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 3.151504E+00 2.051857E+00 -4.923938E-02 @@ -160,6 +310,56 @@ tally 2: 1.674419E-04 -1.682320E-02 8.389254E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 4.536316E+01 4.258781E+02 -6.747275E-01 diff --git a/tests/test_score_total_yn/tallies.xml b/tests/test_score_total_yn/tallies.xml index df6de1c134..ca5dcd2dcb 100644 --- a/tests/test_score_total_yn/tallies.xml +++ b/tests/test_score_total_yn/tallies.xml @@ -9,6 +9,7 @@ total-y4 + U-235 total diff --git a/tests/test_union_energy_grids/geometry.xml b/tests/test_union_energy_grids/geometry.xml new file mode 100644 index 0000000000..612e46132e --- /dev/null +++ b/tests/test_union_energy_grids/geometry.xml @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/tests/test_union_energy_grids/materials.xml b/tests/test_union_energy_grids/materials.xml new file mode 100644 index 0000000000..45ccc65540 --- /dev/null +++ b/tests/test_union_energy_grids/materials.xml @@ -0,0 +1,11 @@ + + + + + + + + + + + diff --git a/tests/test_union_energy_grids/results.py b/tests/test_union_energy_grids/results.py new file mode 100644 index 0000000000..be13ee66f1 --- /dev/null +++ b/tests/test_union_energy_grids/results.py @@ -0,0 +1,25 @@ +#!/usr/bin/env python + +import sys + +# import statepoint +sys.path.insert(0, '../../src/utils') +import statepoint + +# read in statepoint file +if len(sys.argv) > 1: + sp = statepoint.StatePoint(sys.argv[1]) +else: + sp = statepoint.StatePoint('statepoint.10.binary') +sp.read_results() + +# set up output string +outstr = '' + +# write out k-combined +outstr += 'k-combined:\n' +outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1]) + +# write results to file +with open('results_test.dat','w') as fh: + fh.write(outstr) diff --git a/tests/test_union_energy_grids/results_true.dat b/tests/test_union_energy_grids/results_true.dat new file mode 100644 index 0000000000..05cda2e980 --- /dev/null +++ b/tests/test_union_energy_grids/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +3.215828E-01 2.966835E-03 diff --git a/tests/test_union_energy_grids/settings.xml b/tests/test_union_energy_grids/settings.xml new file mode 100644 index 0000000000..1eb22241cc --- /dev/null +++ b/tests/test_union_energy_grids/settings.xml @@ -0,0 +1,18 @@ + + + + union + + + 10 + 5 + 1000 + + + + + -4 -4 -4 4 4 4 + + + + diff --git a/tests/test_union_energy_grids/test_union_energy_grids.py b/tests/test_union_energy_grids/test_union_energy_grids.py new file mode 100644 index 0000000000..6fdbf87459 --- /dev/null +++ b/tests/test_union_energy_grids/test_union_energy_grids.py @@ -0,0 +1,59 @@ +#!/usr/bin/env python + +import os +from subprocess import Popen, STDOUT, PIPE, call +import filecmp +import glob +from optparse import OptionParser + +parser = OptionParser() +parser.add_option('--mpi_exec', dest='mpi_exec', default='') +parser.add_option('--mpi_np', dest='mpi_np', default='3') +parser.add_option('--exe', dest='exe') +(opts, args) = parser.parse_args() +cwd = os.getcwd() + +def test_run(): + if opts.mpi_exec != '': + proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd], + stderr=STDOUT, stdout=PIPE) + else: + proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE) + print(proc.communicate()[0]) + returncode = proc.returncode + assert returncode == 0, 'OpenMC did not exit successfully.' + +def test_created_statepoint(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.' + assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ + 'Statepoint file is not a binary or hdf5 file.' + +def test_results(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + call(['python', 'results.py', statepoint[0]]) + compare = filecmp.cmp('results_test.dat', 'results_true.dat') + if not compare: + os.rename('results_test.dat', 'results_error.dat') + assert compare, 'Results do not agree.' + +def teardown(): + output = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + output.append(os.path.join(cwd, 'results_test.dat')) + for f in output: + if os.path.exists(f): + os.remove(f) + +if __name__ == '__main__': + + # test for openmc executable + if opts.exe is None: + raise Exception('Must specify OpenMC executable from command line with --exe.') + + # run tests + try: + test_run() + test_created_statepoint() + test_results() + finally: + teardown()