From 5fb7d5c2386731ee4e28c1642544cf9515e3f615 Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Sat, 19 Jul 2014 16:09:35 -0600 Subject: [PATCH 01/29] implemented a unionized material energy grid option --- src/ace_header.F90 | 6 +- src/constants.F90 | 7 ++- src/cross_section.F90 | 84 ++++++++++++++++--------- src/energy_grid.F90 | 136 ++++++++++++++++++++++++++++------------ src/initialize.F90 | 4 +- src/input_xml.F90 | 8 ++- src/material_header.F90 | 7 +++ 7 files changed, 175 insertions(+), 77 deletions(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index e5bf1bb3a4..480397a5e6 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -97,7 +97,7 @@ module ace_header ! Energy grid information integer :: n_grid ! # of nuclide grid points - integer, allocatable :: grid_index(:) ! pointers to union grid + integer, allocatable :: glob_grid_index(:) ! pointers to union grid real(8), allocatable :: energy(:) ! energy values corresponding to xs ! Microscopic cross sections @@ -337,8 +337,8 @@ module ace_header integer :: i ! Loop counter - if (allocated(this % grid_index)) & - deallocate(this % grid_index) + if (allocated(this % glob_grid_index)) & + deallocate(this % glob_grid_index) if (allocated(this % energy)) & deallocate(this % total, this % elastic, this % fission, & diff --git a/src/constants.F90 b/src/constants.F90 index 94edd9d2a2..4cedf798c5 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -382,9 +382,10 @@ module constants ! Energy grid methods integer, parameter :: & - GRID_NUCLIDE = 1, & ! non-unionized energy grid - GRID_UNION = 2, & ! union grid with pointers - GRID_LETHARGY = 3 ! lethargy mapping + GRID_NUCLIDE = 1, & ! non-unionized energy grid + GRID_GLOB_UNION = 2, & ! global union grid with pointers + GRID_MAT_UNION = 3, & ! material union grids with pointers + GRID_LETHARGY = 4 ! lethargy mapping ! Running modes integer, parameter :: & diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 7426519cda..8fd2fac4d0 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -37,11 +37,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 @@ -49,8 +49,9 @@ contains mat => materials(p % material) - ! Find energy index on unionized grid - if (grid_method == GRID_UNION) call find_energy_index(p % E) + ! Find energy index on global or material unionized grid + if (grid_method == GRID_GLOB_UNION .or. 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) @@ -92,9 +93,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 ! ======================================================================== @@ -135,34 +136,42 @@ 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 + integer, intent(in) :: i_mat ! index into materials array + integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material real(8), intent(in) :: E ! energy integer :: i_grid ! index on nuclide energy grid real(8) :: f ! interp factor on nuclide energy grid - type(Nuclide), pointer, save :: nuc => null() + type(Nuclide), pointer, save :: nuc => null() + type(Material), pointer, save :: mat => null() !$omp threadprivate(nuc) - ! 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_UNION) - ! If we're using the unionized grid with pointers, finding the index on - ! the nuclide energy grid is as simple as looking up the pointer + ! If we're using a unionized grid with pointers, finding the index on + ! the nuclide energy grid is as simple as looking up the pointer + case (GRID_GLOB_UNION) - i_grid = nuc % grid_index(union_grid_index) + i_grid = nuc % glob_grid_index(union_grid_index) + + case (GRID_MAT_UNION) + + i_grid = mat % nuclide_grid_index(i_nuc_mat, union_grid_index) case (GRID_NUCLIDE) ! If we're not using the unionized grid, we have to do a 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 @@ -465,19 +474,38 @@ contains ! energy !=============================================================================== - subroutine find_energy_index(E) + subroutine find_energy_index(E, i_mat) - real(8), intent(in) :: E ! energy of particle + 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 - ! if particle's 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 < e_grid(1)) then - union_grid_index = 1 - elseif (E > e_grid(n_grid)) then - union_grid_index = n_grid - 1 - else - union_grid_index = binary_search(e_grid, n_grid, E) - end if + select case(grid_method) + case(GRID_GLOB_UNION) + ! 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 < e_grid(1)) then + union_grid_index = 1 + elseif (E > e_grid(n_grid)) then + union_grid_index = n_grid - 1 + else + union_grid_index = binary_search(e_grid, n_grid, E) + end if + + case(GRID_MAT_UNION) + 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 select end subroutine find_energy_index diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 66b95be0de..1c86982fcf 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -8,43 +8,75 @@ module energy_grid contains !=============================================================================== -! UNIONIZED_GRID creates a single unionized energy grid combined from each -! nuclide of each material. 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?) +! 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 :: list => null() - type(Nuclide), pointer :: nuc => null() + type(Nuclide), pointer :: nuc => null() + type(Material), pointer :: mat => null() - message = "Creating unionized energy grid..." + message = "Creating unionized energy grid(s)..." call write_message(5) - ! Add grid points for each nuclide in the problem - do i = 1, n_nuclides_total - nuc => nuclides(i) - call add_grid_points(list, nuc % energy) - end do + select case(grid_method) + case(GRID_GLOB_UNION) + ! Add grid points for each nuclide in the problem + do i = 1, n_nuclides_total + nuc => nuclides(i) + call add_grid_points(list, nuc % energy) + end do - ! Set size of unionized energy grid - n_grid = list % size() + ! Set size of unionized energy grid + n_grid = list % size() - ! create allocated array from linked list - allocate(e_grid(n_grid)) - do i = 1, n_grid - e_grid(i) = list % get_item(i) - end do + ! create allocated array from linked list + allocate(e_grid(n_grid)) + do i = 1, n_grid + e_grid(i) = list % get_item(i) + end do - ! delete linked list and dictionary - call list % clear() - deallocate(list) + ! delete linked list and dictionary + call list % clear() + deallocate(list) - ! Set pointers to unionized energy grid for each nuclide - call grid_pointers() + ! Set pointers to unionized energy grid for each nuclide + call grid_pointers() + + case(GRID_MAT_UNION) + ! 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 select end subroutine unionized_grid @@ -117,34 +149,60 @@ contains !=============================================================================== ! GRID_POINTERS creates an array of pointers (ints) for each nuclide to link -! each point on the nuclide energy grid to one on the unionized energy grid +! 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(Nuclide), pointer :: nuc => null() + type(Material), pointer :: mat => null() - do i = 1, n_nuclides_total - nuc => nuclides(i) - allocate(nuc % grid_index(n_grid)) + select case(grid_method) + case(GRID_GLOB_UNION) + do i = 1, n_nuclides_total + nuc => nuclides(i) + allocate(nuc % glob_grid_index(n_grid)) - index_e = 1 - energy = nuc % energy(index_e) + index_e = 1 + energy = nuc % energy(index_e) - do j = 1, n_grid - union_energy = 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 - nuc % grid_index(j) = index_e - 1 + do j = 1, n_grid + union_energy = 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 + nuc % glob_grid_index(j) = index_e - 1 + end do end do - end do + + case(GRID_MAT_UNION) + 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 select end subroutine grid_pointers diff --git a/src/initialize.F90 b/src/initialize.F90 index c14a8a02a7..bfcd1cc5d7 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -105,7 +105,9 @@ contains call time_read_xs % stop() ! Construct unionized energy grid from cross-sections - if (grid_method == GRID_UNION) then + if (grid_method == GRID_NUCLIDE) then + continue + else call time_unionize % start() call unionized_grid() call time_unionize % stop() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 46461b78da..b6523f6eab 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -213,13 +213,15 @@ contains if (check_for_node(doc, "energy_grid")) then call get_node_value(doc, "energy_grid", temp_str) else - temp_str = 'union' + temp_str = 'global-union' end if select case (trim(temp_str)) case ('nuclide') grid_method = GRID_NUCLIDE - case ('union') - grid_method = GRID_UNION + case ('global-union') + grid_method = GRID_GLOB_UNION + case ('material-union') + grid_method = GRID_MAT_UNION case ('lethargy') message = "Lethargy mapped energy grid not yet supported." call fatal_error() diff --git a/src/material_header.F90 b/src/material_header.F90 index cbfd917498..6b6bb593eb 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 From 69f378518c55cdc81744690322cc96ab205f1dbd Mon Sep 17 00:00:00 2001 From: walshjon Date: Thu, 13 Nov 2014 20:10:57 -0800 Subject: [PATCH 02/29] remove unnecessary critical blocks to speed up threading --- src/cross_section.F90 | 3 ++- src/energy_grid.F90 | 7 ++++--- src/initialize.F90 | 2 ++ 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index e83e1303a4..370c5c3630 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -149,7 +149,7 @@ contains real(8) :: f ! interp factor on nuclide energy grid type(Nuclide), pointer, save :: nuc => null() type(Material), pointer, save :: mat => null() -!$omp threadprivate(nuc) +!$omp threadprivate(nuc, mat) ! Set pointer to nuclide and material nuc => nuclides(i_nuclide) @@ -515,6 +515,7 @@ contains 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) select case(grid_method) case(GRID_GLOB_UNION) diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 7ea9d34ff4..b168d17ff1 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -22,9 +22,10 @@ contains integer :: i ! index in nuclides array integer :: j ! index in materials array - type(ListReal), pointer :: list => null() - type(Nuclide), pointer :: nuc => null() - type(Material), pointer :: mat => null() + 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) diff --git a/src/initialize.F90 b/src/initialize.F90 index 34272fe882..f15541df0f 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -113,7 +113,9 @@ contains continue else call time_unionize % start() +!$omp critical call unionized_grid() +!$omp end critical call time_unionize % stop() end if From 3858f412374c1200199262633af525179c599eee Mon Sep 17 00:00:00 2001 From: walshjon Date: Thu, 19 Mar 2015 18:19:37 -0700 Subject: [PATCH 03/29] cleaned up merge of mit-crpg/openmc/develop --- src/ace_header.F90 | 3 -- src/constants.F90 | 2 +- src/cross_section.F90 | 3 ++ src/energy_grid.F90 | 71 ++++++++++++++++++++++++++++++++++++++++++- src/global.F90 | 1 + src/initialize.F90 | 4 +-- 6 files changed, 77 insertions(+), 7 deletions(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 655b1480ec..7753f4734a 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -372,9 +372,6 @@ module ace_header integer :: i ! Loop counter - if (allocated(this % glob_grid_index)) & - deallocate(this % glob_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 bd8d805783..85976c00ac 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -368,7 +368,7 @@ module constants ! Energy grid methods integer, parameter :: & - GRID_NUCLIDE = 1, & ! non-unionized energy grid + GRID_NUCLIDE = 1, & ! unique energy grid for each nuclide GRID_MAT_UNION = 2, & ! material union grids with pointers GRID_LOGARITHM = 3 ! lethargy mapping diff --git a/src/cross_section.F90 b/src/cross_section.F90 index c667dd2a52..8a6981d223 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 !=============================================================================== diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 5ac8f561ed..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 @@ -56,7 +59,6 @@ contains call grid_pointers() end subroutine unionized_grid ->>>>>>> mit-crpg/openmc/develop !=============================================================================== ! LOGARITHMIC_GRID determines a logarithmic mapping for energies to bounding @@ -107,6 +109,73 @@ 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 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 8c5079c789..29ed52aa0b 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,& @@ -121,7 +121,7 @@ contains 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() From 1e2f803776e525ed5f20821b417262603754f2ad Mon Sep 17 00:00:00 2001 From: walshjon Date: Thu, 26 Mar 2015 13:19:48 -0700 Subject: [PATCH 04/29] add test, address review comments --- src/cross_section.F90 | 2 +- src/initialize.F90 | 2 - src/input_xml.F90 | 5 +- tests/test_union_energy_grids/geometry.xml | 8 +++ tests/test_union_energy_grids/materials.xml | 11 ++++ tests/test_union_energy_grids/results.py | 25 ++++++++ .../test_union_energy_grids/results_true.dat | 2 + tests/test_union_energy_grids/settings.xml | 18 ++++++ .../test_union_energy_grids.py | 59 +++++++++++++++++++ 9 files changed, 128 insertions(+), 4 deletions(-) create mode 100644 tests/test_union_energy_grids/geometry.xml create mode 100644 tests/test_union_energy_grids/materials.xml create mode 100644 tests/test_union_energy_grids/results.py create mode 100644 tests/test_union_energy_grids/results_true.dat create mode 100644 tests/test_union_energy_grids/settings.xml create mode 100644 tests/test_union_energy_grids/test_union_energy_grids.py diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 8a6981d223..dce4d75f32 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -53,7 +53,7 @@ contains ! Find energy index on global or material unionized grid if (grid_method == GRID_MAT_UNION) & - & call find_energy_index(p % E, p % material) + call find_energy_index(p % E, p % material) ! Determine if this material has S(a,b) tables check_sab = (mat % n_sab > 0) diff --git a/src/initialize.F90 b/src/initialize.F90 index 29ed52aa0b..51dbf5ed52 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -115,9 +115,7 @@ contains continue case (GRID_MAT_UNION) call time_unionize % start() -!$omp critical call unionized_grid() -!$omp end critical call time_unionize % stop() case (GRID_LOGARITHM) call logarithmic_grid() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 0ae299aa4e..dd28cd9610 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 ('material-union') + 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/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() From 66e669f07671c5ab0461c7a09f4eacb3bf2c8dbc Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 6 Apr 2015 21:29:02 -0400 Subject: [PATCH 05/29] Fixed tracklength total yn score --- src/tally.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/tally.F90 b/src/tally.F90 index c96f1ecd6e..39067a9c07 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -782,6 +782,9 @@ contains 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 From 7ac15772e3b507ed3344c7f3e64905978c1c4031 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 7 Apr 2015 12:40:57 -0400 Subject: [PATCH 06/29] Added general tally scoring subroutine --- src/tally.F90 | 1296 +++++++++++++++++++++++++++++++++---------------- 1 file changed, 873 insertions(+), 423 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index c96f1ecd6e..a3d647573c 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -26,6 +26,451 @@ module tally contains + 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 + real(8), intent(in) :: flux + real(8), intent(in) :: atom_density ! atom/b-cm + + integer :: j ! 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) :: score_ ! analog tally score + real(8) :: macro_total ! material macro total xs + real(8) :: macro_scatt ! material macro scatt xs + type(Material), pointer, save :: mat => null() + type(Reaction), pointer, save :: rxn => null() +!$omp threadprivate(t, mat, rxn) + + 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 = start_index + j + + + 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 * atom_density * 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) * atom_density * 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 + j = j + t % moment_order(j) + 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 + j = j + t % moment_order(j) + 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 * & + µ_xs(p % event_nuclide) % kappa_fission / & + µ_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 * & + µ_xs(p % event_nuclide) % kappa_fission / & + µ_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 + + + ! Add score to tally. + 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(j) == 1) then + score = score * p % mu ! avoid function call overhead + else + score = score * calc_pn(t % moment_order(j), p % mu) + endif +!$omp atomic + t % results(score_index, filter_index) % value = & + &t % results(score_index, filter_index) % value + score + + + case(SCORE_FLUX_YN, SCORE_TOTAL_YN, 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(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_pn(n, p % mu) * calc_rn(n, p % last_uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 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(j) + ! determine scoring bin index + score_index = score_index + 1 + ! get the score and tally it + score_ = score * calc_pn(n, p % 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) + + + 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 @@ -122,429 +567,434 @@ contains ! Now compare the value against that of the colliding nuclide. if (i_nuclide /= p % event_nuclide .and. i_nuclide /= -1) cycle 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 +! subroutine score_general(p, t, start_index, filter_index, i_nuclide, & +! &atom_density, flux) + + call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & + &i_nuclide, ZERO, ZERO) + +! ! 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 end do NUCLIDE_LOOP From ac923dd8e154f24ae56a0cd4387e3dfbef4940df Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 7 Apr 2015 16:22:29 -0400 Subject: [PATCH 07/29] Added tracklength tallies to the generalized score --- src/tally.F90 | 609 ++++++++++++++++++++++++++------------------------ 1 file changed, 320 insertions(+), 289 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index a3d647573c..f64fda7c18 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -52,6 +52,7 @@ contains 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(t, mat, rxn) @@ -106,7 +107,7 @@ contains if (i_nuclide > 0) then score = micro_xs(i_nuclide) % total * atom_density * flux else - score = material_xs % total * atom_density * flux + score = material_xs % total * flux end if end if @@ -126,8 +127,7 @@ contains score = (micro_xs(i_nuclide) % total & &- micro_xs(i_nuclide) % absorption) * atom_density * flux else - score = (material_xs % total & - &- material_xs % absorption) * atom_density * flux + score = (material_xs % total - material_xs % absorption) * flux end if end if @@ -423,7 +423,7 @@ contains &t % results(score_index, filter_index) % value + score - case(SCORE_FLUX_YN, SCORE_TOTAL_YN, SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN) + 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 @@ -444,6 +444,32 @@ contains j = j + (t % moment_order(j) + 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(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, uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 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 @@ -1179,293 +1205,298 @@ contains atom_density = ZERO end if end if +! subroutine score_general(p, t, start_index, filter_index, i_nuclide, & +! &atom_density, flux) - ! 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) + call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & + &i_nuclide, atom_density, flux) - ! 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 +! ! 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 end do NUCLIDE_BIN_LOOP end if From 347610a15329d9f24a62301e99363c3c048a5a6c Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 7 Apr 2015 17:14:09 -0400 Subject: [PATCH 08/29] Added tracklength mesh to general scoring --- src/tally.F90 | 926 +------------------------------------------------- 1 file changed, 11 insertions(+), 915 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index f64fda7c18..fb0b95a5b8 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -26,6 +26,12 @@ 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 @@ -509,22 +515,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) @@ -593,435 +590,11 @@ contains ! Now compare the value against that of the colliding nuclide. if (i_nuclide /= p % event_nuclide .and. i_nuclide /= -1) cycle end if -! subroutine score_general(p, t, start_index, filter_index, i_nuclide, & -! &atom_density, flux) + ! Determine score for each bin call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & &i_nuclide, ZERO, ZERO) -! ! 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 - end do NUCLIDE_LOOP ! If the user has specified that we can assume all tallies are spatially @@ -1116,26 +689,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 @@ -1205,299 +766,11 @@ contains atom_density = ZERO end if end if -! subroutine score_general(p, t, start_index, filter_index, i_nuclide, & -! &atom_density, flux) + ! Determine score for each bin call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & &i_nuclide, atom_density, flux) -! ! 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 - end do NUCLIDE_BIN_LOOP end if @@ -1865,21 +1138,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 @@ -2093,179 +1360,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 From f24e14fbd531e6b7dfa38da1380f929deab790bf Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 7 Apr 2015 17:37:34 -0400 Subject: [PATCH 09/29] Added general scoring to score_all_nuclides --- src/tally.F90 | 298 ++------------------------------------------------ 1 file changed, 9 insertions(+), 289 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index fb0b95a5b8..2e59301e1c 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -801,22 +801,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) @@ -835,290 +824,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 From d9e0f6cee81b084f33ccf7fc664c67a5dc8b30a3 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 9 Apr 2015 12:31:58 -0400 Subject: [PATCH 10/29] Added nuclides to test_score_total_yn --- tests/test_score_total_yn/results_true.dat | 200 +++++++++++++++++++++ tests/test_score_total_yn/tallies.xml | 1 + 2 files changed, 201 insertions(+) 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 From f63c3461d7a6e1f428c8cf1bc23a1a5fcfd16bd6 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 9 Apr 2015 22:21:58 -0400 Subject: [PATCH 11/29] Imporved score_general readability --- src/tally.F90 | 64 +++++++++++++++++++++++++++++---------------------- 1 file changed, 36 insertions(+), 28 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 2e59301e1c..faea7e2ebd 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -38,11 +38,11 @@ contains type(TallyObject), pointer, intent(inout) :: t integer, intent(in) :: start_index integer, intent(in) :: i_nuclide - integer, intent(in) :: filter_index - real(8), intent(in) :: flux + integer, intent(in) :: filter_index ! for % results + real(8), intent(in) :: flux ! flux estimate real(8), intent(in) :: atom_density ! atom/b-cm - integer :: j ! loop index for scoring bins + 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 @@ -63,15 +63,18 @@ contains type(Reaction), pointer, save :: rxn => null() !$omp threadprivate(t, mat, rxn) - j = 0 + i = 0 SCORE_LOOP: do q = 1, t % n_user_score_bins - j = j + 1 + i = i + 1 + ! determine what type of score bin - score_bin = t % score_bins(j) + score_bin = t % score_bins(i) ! determine scoring bin index - score_index = start_index + j + score_index = start_index + i + !######################################################################### + ! Determine appropirate scoring value. select case(score_bin) @@ -142,7 +145,7 @@ contains ! Only analog estimators are available. ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) then - j = j + t % moment_order(j) + i = i + t % moment_order(i) cycle SCORE_LOOP end if ! Since only scattering events make it here, again we can use @@ -165,7 +168,7 @@ contains ! Only analog estimators are available. ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) then - j = j + t % moment_order(j) + i = i + t % moment_order(i) cycle SCORE_LOOP end if ! For scattering production, we need to use the post-collision @@ -185,7 +188,7 @@ contains ! 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) + score = (macro_total - p % mu * macro_scatt) * (ONE / macro_scatt) case (SCORE_N_1N) @@ -411,18 +414,19 @@ contains end select + !######################################################################### + ! Expand score if necessary and add to tally results. - ! Add score to tally. 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(j) == 1) then + if (t % moment_order(i) == 1) then score = score * p % mu ! avoid function call overhead else - score = score * calc_pn(t % moment_order(j), p % mu) + score = score * calc_pn(t % moment_order(i), p % mu) endif !$omp atomic t % results(score_index, filter_index) % value = & @@ -434,7 +438,7 @@ contains 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) + 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]) @@ -442,12 +446,14 @@ contains ! 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) + 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 - j = j + (t % moment_order(j) + 1)**2 - 1 + i = i + (t % moment_order(i) + 1)**2 - 1 case(SCORE_FLUX_YN, SCORE_TOTAL_YN) @@ -460,7 +466,7 @@ contains end if ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, t % moment_order(j) + 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]) @@ -468,29 +474,31 @@ contains ! 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) + 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 - j = j + (t % moment_order(j) + 1)**2 - 1 + 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(j) + do n = 0, t % moment_order(i) ! determine scoring bin index score_index = score_index + 1 - ! get the score and tally it - score_ = score * calc_pn(n, p % mu) + ! get the score and tally it !$omp atomic t % results(score_index, filter_index) % value = & - &t % results(score_index, filter_index) % value + score_ + & t % results(score_index, filter_index) % value & + & + score * calc_pn(n, p % mu) end do - j = j + t % moment_order(j) + i = i + t % moment_order(i) case default From 37e00a484f292af0223b686f7af364d466f9408e Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 13 Apr 2015 16:42:05 -0400 Subject: [PATCH 12/29] Fixed OpenMP for score_general --- src/tally.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index faea7e2ebd..fd2157e810 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -55,13 +55,12 @@ contains real(8) :: atom_density_ ! atom/b-cm real(8) :: f ! interpolation factor real(8) :: score ! analog tally score - 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(t, mat, rxn) +!$omp threadprivate(mat, rxn) i = 0 SCORE_LOOP: do q = 1, t % n_user_score_bins From d2d386f5a43328c72b3cb24b365048450780605e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 14 Apr 2015 10:14:30 -0500 Subject: [PATCH 13/29] Better recognition of MPI wrappers in CMakeLists.txt --- src/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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) From 2abe82c696b9ffb8610500cbbf0ed1c7d6899d5f Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 15 Apr 2015 11:31:19 -0400 Subject: [PATCH 14/29] Removed ampersands for #370 --- src/tally.F90 | 54 +++++++++++++++++++++++++-------------------------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index fd2157e810..ded646d4da 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -33,7 +33,7 @@ contains !=============================================================================== subroutine score_general(p, t, start_index, filter_index, i_nuclide, & - &atom_density, flux) + atom_density, flux) type(Particle), intent(in) :: p type(TallyObject), pointer, intent(inout) :: t integer, intent(in) :: start_index @@ -133,7 +133,7 @@ contains ! 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 + - micro_xs(i_nuclide) % absorption) * atom_density * flux else score = (material_xs % total - material_xs % absorption) * flux end if @@ -231,7 +231,7 @@ contains ! 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 + / micro_xs(p % event_nuclide) % absorption else score = ZERO end if @@ -242,7 +242,7 @@ contains ! 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 + / micro_xs(p % event_nuclide) % absorption end if else if (t % estimator == ESTIMATOR_TRACKLENGTH) then @@ -273,7 +273,7 @@ contains ! 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 + nu_fission / micro_xs(p % event_nuclide) % absorption else score = ZERO end if @@ -305,8 +305,8 @@ contains ! fission scale by kappa-fission if (micro_xs(p % event_nuclide) % absorption > ZERO) then score = p % absorb_wgt * & - µ_xs(p % event_nuclide) % kappa_fission / & - µ_xs(p % event_nuclide) % absorption + micro_xs(p % event_nuclide) % kappa_fission / & + micro_xs(p % event_nuclide) % absorption else score = ZERO end if @@ -317,8 +317,8 @@ contains ! particle's weight entering the collision as the estimate for ! the fission energy production rate score = p % last_wgt * & - µ_xs(p % event_nuclide) % kappa_fission / & - µ_xs(p % event_nuclide) % absorption + micro_xs(p % event_nuclide) % kappa_fission / & + micro_xs(p % event_nuclide) % absorption end if else if (t % estimator == ESTIMATOR_TRACKLENGTH) then @@ -406,7 +406,7 @@ contains else call fatal_error("Invalid score type on tally " & - &// to_str(t % id) // ".") + // to_str(t % id) // ".") end if end if @@ -429,7 +429,7 @@ contains endif !$omp atomic t % results(score_index, filter_index) % value = & - &t % results(score_index, filter_index) % value + score + t % results(score_index, filter_index) % value + score case(SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN) @@ -446,10 +446,10 @@ contains ! 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) + % 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 @@ -474,10 +474,10 @@ contains ! 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) + % 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 @@ -494,8 +494,8 @@ contains ! 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) + t % results(score_index, filter_index) % value & + + score * calc_pn(n, p % mu) end do i = i + t % moment_order(i) @@ -503,7 +503,7 @@ contains case default !$omp atomic t % results(score_index, filter_index) % value = & - &t % results(score_index, filter_index) % value + score + t % results(score_index, filter_index) % value + score end select @@ -600,7 +600,7 @@ contains ! Determine score for each bin call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & - &i_nuclide, ZERO, ZERO) + i_nuclide, ZERO, ZERO) end do NUCLIDE_LOOP @@ -776,7 +776,7 @@ contains ! Determine score for each bin call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & - &i_nuclide, atom_density, flux) + i_nuclide, atom_density, flux) end do NUCLIDE_BIN_LOOP end if @@ -833,7 +833,7 @@ contains ! Determine score for each bin call score_general(p, t, (i_nuclide-1)*t % n_score_bins, filter_index, & - &i_nuclide, atom_density, flux) + i_nuclide, atom_density, flux) end do NUCLIDE_LOOP @@ -845,7 +845,7 @@ contains ! Determine score for each bin call score_general(p, t, n_nuclides_total*t % n_score_bins, filter_index, & - &i_nuclide, atom_density, flux) + i_nuclide, atom_density, flux) end subroutine score_all_nuclides @@ -1088,7 +1088,7 @@ contains ! Determine score for each bin call score_general(p, t, (b-1)*t % n_score_bins, filter_index, & - &i_nuclide, atom_density, flux) + i_nuclide, atom_density, flux) end do NUCLIDE_BIN_LOOP end if From 500fea92cc082c7ddff31a256aa2c9a59e931f0b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 15 Apr 2015 15:00:14 -0500 Subject: [PATCH 15/29] Make sure MPI_DIR and PETSC_DIR are passed when running CMake from run_tests.py --- tests/run_tests.py | 4 ++++ 1 file changed, 4 insertions(+) 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) From 49b09c00dae175e57a701c045a577470d69c816b Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Thu, 16 Apr 2015 01:42:58 -0400 Subject: [PATCH 16/29] update schema for settings.xml --- src/relaxng/settings.rnc | 2 +- src/relaxng/settings.rng | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) 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 From 057d65c2b1a00cdfc0f298e3ceea7c47904c6ba5 Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Thu, 16 Apr 2015 01:54:18 -0400 Subject: [PATCH 17/29] update documentation to mention energy grids unionized by material --- docs/source/usersguide/input.rst | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 79d66d2869..2ca93e3288 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 From f0a3905b1bfabd522b8e63dd1f9e97516d01dd4e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Apr 2015 10:13:57 -0500 Subject: [PATCH 18/29] Don't show k-effective in results when run in fixed source mode --- src/output.F90 | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) 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,*) From 120b283e577fa1da59f1f8a3f13469f60ea9db88 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Apr 2015 10:26:34 -0500 Subject: [PATCH 19/29] Update acceptable options for energy_grid --- src/utils/openmc/settings.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/utils/openmc/settings.py b/src/utils/openmc/settings.py index b5e195e3f7..0e314545a2 100644 --- a/src/utils/openmc/settings.py +++ b/src/utils/openmc/settings.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - import warnings from openmc.checkvalue import * @@ -437,9 +435,9 @@ class SettingsFile(object): def set_energy_grid(self, energy_grid): - if not energy_grid in ['union', 'nuclide']: + if not energy_grid in ['nuclide', 'logarithm', 'material-union']: msg = 'Unable to set energy grid to {0} which is neither ' \ - 'union nor nuclide'.format(energy_grid) + 'nuclide, logarithm, nor material-union'.format(energy_grid) raise ValueError(msg) self._energy_grid = energy_grid @@ -1206,7 +1204,7 @@ class SettingsFile(object): subelement.text = 'true' else: subelement.text = 'false' - + subelement = ET.SubElement(element, "count_interactions") if self._dd_count_interactions: subelement.text = 'true' From 2d92f81ed9d634621260c4726d1f52a03c0cd7cf Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Apr 2015 12:27:27 -0500 Subject: [PATCH 20/29] Ignore example inputs generated by Python API --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index e758f64e43..a5cb06dc4f 100644 --- a/.gitignore +++ b/.gitignore @@ -11,6 +11,9 @@ # OpenMC executable src/openmc +# Inputs generated from Python API +examples/python/**/*.xml + # emacs backups *~ From 764ab84f0136484adbc94b75bc2e29373991e5e7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Apr 2015 12:28:04 -0500 Subject: [PATCH 21/29] Remove #! from Python files in Python API since they don't have __main__ --- src/utils/openmc/__init__.py | 2 -- src/utils/openmc/checkvalue.py | 3 --- src/utils/openmc/clean_xml.py | 4 +--- src/utils/openmc/cmfd.py | 4 +--- src/utils/openmc/constants.py | 2 -- src/utils/openmc/element.py | 2 -- src/utils/openmc/executor.py | 4 +--- src/utils/openmc/geometry.py | 4 +--- src/utils/openmc/material.py | 4 +--- src/utils/openmc/nuclide.py | 4 +--- src/utils/openmc/opencg_compatible.py | 2 -- src/utils/openmc/plots.py | 2 -- src/utils/openmc/statepoint.py | 4 +--- src/utils/openmc/summary.py | 4 +--- src/utils/openmc/surface.py | 4 +--- src/utils/openmc/tallies.py | 2 -- src/utils/openmc/universe.py | 2 -- 17 files changed, 9 insertions(+), 44 deletions(-) diff --git a/src/utils/openmc/__init__.py b/src/utils/openmc/__init__.py index 2f3a9b6de4..42ee1eebfa 100644 --- a/src/utils/openmc/__init__.py +++ b/src/utils/openmc/__init__.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - from openmc.element import * from openmc.geometry import * from openmc.nuclide import * diff --git a/src/utils/openmc/checkvalue.py b/src/utils/openmc/checkvalue.py index 679b5c2e42..503e7e24d5 100644 --- a/src/utils/openmc/checkvalue.py +++ b/src/utils/openmc/checkvalue.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - import numpy as np @@ -13,4 +11,3 @@ def is_float(val): def is_string(val): return isinstance(val, (str, np.str)) - diff --git a/src/utils/openmc/clean_xml.py b/src/utils/openmc/clean_xml.py index cac5833dc3..d67c01d549 100644 --- a/src/utils/openmc/clean_xml.py +++ b/src/utils/openmc/clean_xml.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - def sort_xml_elements(tree): # Retrieve all children of the root XML node in the tree @@ -92,4 +90,4 @@ def clean_xml_indentation(element, level=0): else: if level and (not element.tail or not element.tail.strip()): - element.tail = i \ No newline at end of file + element.tail = i diff --git a/src/utils/openmc/cmfd.py b/src/utils/openmc/cmfd.py index 4b7802ba3d..443867284a 100644 --- a/src/utils/openmc/cmfd.py +++ b/src/utils/openmc/cmfd.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - from openmc.checkvalue import * from openmc.clean_xml import * from xml.etree import ElementTree as ET @@ -583,4 +581,4 @@ class CMFDFile(object): # Write the XML Tree to the cmfd.xml file tree = ET.ElementTree(self._cmfd_file) tree.write("cmfd.xml", xml_declaration=True, - encoding='utf-8', method="xml") \ No newline at end of file + encoding='utf-8', method="xml") diff --git a/src/utils/openmc/constants.py b/src/utils/openmc/constants.py index ef55d2f8f2..7c4a7aed41 100644 --- a/src/utils/openmc/constants.py +++ b/src/utils/openmc/constants.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - """Dictionaries of integer-to-string mappings from openmc/src/constants.F90""" SURFACE_TYPES = {1: 'x-plane', diff --git a/src/utils/openmc/element.py b/src/utils/openmc/element.py index 745e2ad6bc..4a6a189dc5 100644 --- a/src/utils/openmc/element.py +++ b/src/utils/openmc/element.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - from openmc.checkvalue import * class Element(object): diff --git a/src/utils/openmc/executor.py b/src/utils/openmc/executor.py index cb09ea8515..f2d26de9b6 100644 --- a/src/utils/openmc/executor.py +++ b/src/utils/openmc/executor.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - from openmc.checkvalue import * import subprocess import os @@ -72,4 +70,4 @@ class Executor(object): cwd=self._working_directory) else: subprocess.check_call(command, shell=True, stdout=FNULL, - cwd=self._working_directory) \ No newline at end of file + cwd=self._working_directory) diff --git a/src/utils/openmc/geometry.py b/src/utils/openmc/geometry.py index 59fb158ee2..d30da1db51 100644 --- a/src/utils/openmc/geometry.py +++ b/src/utils/openmc/geometry.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - import openmc from openmc.clean_xml import * from xml.etree import ElementTree as ET @@ -159,4 +157,4 @@ class GeometryFile(object): # Write the XML Tree to the materials.xml file tree = ET.ElementTree(self._geometry_file) tree.write("geometry.xml", xml_declaration=True, - encoding='utf-8', method="xml") \ No newline at end of file + encoding='utf-8', method="xml") diff --git a/src/utils/openmc/material.py b/src/utils/openmc/material.py index 914cb154a1..2f501d0763 100644 --- a/src/utils/openmc/material.py +++ b/src/utils/openmc/material.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - import warnings import openmc @@ -365,7 +363,7 @@ class Material(object): else: subelement = ET.SubElement(element, "compositions") - + comps = [] allnucs = self._nuclides.values() + self._elements.values() dist_per_type = allnucs[0][2] diff --git a/src/utils/openmc/nuclide.py b/src/utils/openmc/nuclide.py index 8603c833b8..c351e18081 100644 --- a/src/utils/openmc/nuclide.py +++ b/src/utils/openmc/nuclide.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - from openmc.checkvalue import * @@ -80,4 +78,4 @@ class Nuclide(object): string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs) if self._zaid is not None: string += '{0: <16}{1}{2}\n'.format('\tZAID', '=\t', self._zaid) - return string \ No newline at end of file + return string diff --git a/src/utils/openmc/opencg_compatible.py b/src/utils/openmc/opencg_compatible.py index 2b2323a849..3a45f45005 100644 --- a/src/utils/openmc/opencg_compatible.py +++ b/src/utils/openmc/opencg_compatible.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - import openmc import opencg import copy diff --git a/src/utils/openmc/plots.py b/src/utils/openmc/plots.py index 551b2e76f4..dcb74c8adc 100644 --- a/src/utils/openmc/plots.py +++ b/src/utils/openmc/plots.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - from openmc.checkvalue import * from openmc.clean_xml import * from xml.etree import ElementTree as ET diff --git a/src/utils/openmc/statepoint.py b/src/utils/openmc/statepoint.py index 7c4f90ae38..dad2f4d268 100644 --- a/src/utils/openmc/statepoint.py +++ b/src/utils/openmc/statepoint.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - import struct, copy import numpy as np import scipy.stats @@ -383,7 +381,7 @@ class StatePoint(object): # Extract the moment order string for each score for k in range(len(scores)): - moment = self._get_string(8, + moment = self._get_string(8, path='{0}order{1}'.format(subbase, k+1)) moment = moment.lstrip('[\'') moment = moment.rstrip('\']') diff --git a/src/utils/openmc/summary.py b/src/utils/openmc/summary.py index 81e0b3bf5d..14c06f0460 100644 --- a/src/utils/openmc/summary.py +++ b/src/utils/openmc/summary.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - import openmc from openmc.opencg_compatible import get_opencg_geometry import numpy as np @@ -354,7 +352,7 @@ class Summary(object): if lattice_type == 'rectangular': dimension = self._f['geometry/lattices'][key]['n_cells'][...] lower_left = \ - self._f['geometry/lattices'][key]['lower_left'][...] + self._f['geometry/lattices'][key]['lower_left'][...] pitch = self._f['geometry/lattices'][key]['pitch'][...] outer = self._f['geometry/lattices'][key]['outer'][0] diff --git a/src/utils/openmc/surface.py b/src/utils/openmc/surface.py index 269bfd1c70..087d3981b3 100644 --- a/src/utils/openmc/surface.py +++ b/src/utils/openmc/surface.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - from openmc.checkvalue import * from openmc.constants import BC_TYPES from xml.etree import ElementTree as ET @@ -592,4 +590,4 @@ class ZCone(Cone): # Initialize ZCone class attributes super(ZCone, self).__init__(surface_id, bc_type, x0, y0, z0, R2, name=name) - self._type = 'z-cone' \ No newline at end of file + self._type = 'z-cone' diff --git a/src/utils/openmc/tallies.py b/src/utils/openmc/tallies.py index 063198ea1d..36880236cf 100644 --- a/src/utils/openmc/tallies.py +++ b/src/utils/openmc/tallies.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - from openmc import Nuclide from openmc.clean_xml import * from openmc.checkvalue import * diff --git a/src/utils/openmc/universe.py b/src/utils/openmc/universe.py index d48c6e9ca4..4c5b688b0a 100644 --- a/src/utils/openmc/universe.py +++ b/src/utils/openmc/universe.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - import openmc from openmc.checkvalue import * from xml.etree import ElementTree as ET From eb38ee3b885aa44e3ada82e734c20138c375f4f0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Apr 2015 12:31:26 -0500 Subject: [PATCH 22/29] Use literal empty dict and list instead of dict() and list() --- src/utils/openmc/clean_xml.py | 6 +++--- src/utils/openmc/element.py | 2 +- src/utils/openmc/geometry.py | 4 ++-- src/utils/openmc/material.py | 18 ++++++++-------- src/utils/openmc/nuclide.py | 2 +- src/utils/openmc/opencg_compatible.py | 22 ++++++++++---------- src/utils/openmc/plots.py | 2 +- src/utils/openmc/statepoint.py | 22 ++++++++++---------- src/utils/openmc/summary.py | 20 +++++++++--------- src/utils/openmc/surface.py | 4 ++-- src/utils/openmc/tallies.py | 30 +++++++++++++-------------- src/utils/openmc/universe.py | 24 ++++++++++----------- 12 files changed, 78 insertions(+), 78 deletions(-) diff --git a/src/utils/openmc/clean_xml.py b/src/utils/openmc/clean_xml.py index d67c01d549..437c2df15c 100644 --- a/src/utils/openmc/clean_xml.py +++ b/src/utils/openmc/clean_xml.py @@ -4,7 +4,7 @@ def sort_xml_elements(tree): elements = tree.getchildren() # Initialize empty lists for the sorted and comment elements - sorted_elements = list() + sorted_elements = [] # Initialize an empty set of tags (e.g., Surface, Cell, and Lattice) tags = set() @@ -14,7 +14,7 @@ def sort_xml_elements(tree): tags.add(element.tag) # Initialize an empty list for the comment elements - comment_elements = list() + comment_elements = [] # Find the comment elements and record their ordering within the # tree using a precedence with respect to the subsequent nodes @@ -38,7 +38,7 @@ def sort_xml_elements(tree): continue # Initialize an empty list of tuples to sort (id, element) - tagged_data = list() + tagged_data = [] # Retrieve the IDs for each of the elements for element in tagged_elements: diff --git a/src/utils/openmc/element.py b/src/utils/openmc/element.py index 4a6a189dc5..06b323b12c 100644 --- a/src/utils/openmc/element.py +++ b/src/utils/openmc/element.py @@ -34,7 +34,7 @@ class Element(object): def __hash__(self): - hashable = list() + hashable = [] hashable.append(self._name) hashable.append(self._xs) return hash(tuple(hashable)) diff --git a/src/utils/openmc/geometry.py b/src/utils/openmc/geometry.py index d30da1db51..17d0f5fa87 100644 --- a/src/utils/openmc/geometry.py +++ b/src/utils/openmc/geometry.py @@ -15,7 +15,7 @@ class Geometry(object): # Initialize Geometry class attributes self._root_universe = None - self._offsets = dict() + self._offsets = {} def get_offset(self, path, filter_offset): @@ -61,7 +61,7 @@ class Geometry(object): def get_all_nuclides(self): - nuclides = dict() + nuclides = {} materials = self.get_all_materials() for material in materials: diff --git a/src/utils/openmc/material.py b/src/utils/openmc/material.py index 2f501d0763..0d6e39b992 100644 --- a/src/utils/openmc/material.py +++ b/src/utils/openmc/material.py @@ -10,7 +10,7 @@ import numpy as np # A list of all IDs for all Materials created -MATERIAL_IDS = list() +MATERIAL_IDS = [] # A static variable for auto-generated Material IDs AUTO_MATERIAL_ID = 10000 @@ -18,7 +18,7 @@ AUTO_MATERIAL_ID = 10000 def reset_auto_material_id(): global AUTO_MATERIAL_ID, MATERIAL_IDS AUTO_MATERIAL_ID = 10000 - MATERIAL_IDS = list() + MATERIAL_IDS = [] # Units for density supported by OpenMC @@ -47,15 +47,15 @@ class Material(object): # A dictionary of Nuclides # Keys - Nuclide names # Values - tuple (nuclide, percent, percent type) - self._nuclides = dict() + self._nuclides = {} # A dictionary of Elements # Keys - Element names # Values - tuple (element, percent, percent type) - self._elements = dict() + self._elements = {} # If specified, a list of tuples of (table name, xs identifier) - self._sab = list() + self._sab = [] # If true, the material will be initialized as distributed self._convert_to_distrib_comps = False @@ -241,7 +241,7 @@ class Material(object): def get_all_nuclides(self): - nuclides = dict() + nuclides = {} for nuclide_name, nuclide_tuple in self._nuclides.items(): nuclide = nuclide_tuple[0] @@ -318,7 +318,7 @@ class Material(object): def get_nuclides_xml(self, nuclides, distrib=False): - xml_elements = list() + xml_elements = [] for nuclide in nuclides.values(): xml_elements.append(self.get_nuclide_xml(nuclide, distrib)) @@ -328,7 +328,7 @@ class Material(object): def get_elements_xml(self, elements, distrib=False): - xml_elements = list() + xml_elements = [] for element in elements.values(): xml_elements.append(self.get_element_xml(element, distrib)) @@ -414,7 +414,7 @@ class MaterialsFile(object): def __init__(self): # Initialize MaterialsFile class attributes - self._materials = list() + self._materials = [] self._default_xs = None self._materials_file = ET.Element("materials") diff --git a/src/utils/openmc/nuclide.py b/src/utils/openmc/nuclide.py index c351e18081..e100d6d74e 100644 --- a/src/utils/openmc/nuclide.py +++ b/src/utils/openmc/nuclide.py @@ -36,7 +36,7 @@ class Nuclide(object): def __hash__(self): - hashable = list() + hashable = [] hashable.append(self._name) hashable.append(self._xs) return hash(tuple(hashable)) diff --git a/src/utils/openmc/opencg_compatible.py b/src/utils/openmc/opencg_compatible.py index 3a45f45005..f80eb2c2c7 100644 --- a/src/utils/openmc/opencg_compatible.py +++ b/src/utils/openmc/opencg_compatible.py @@ -7,52 +7,52 @@ import numpy as np # A dictionary of all OpenMC Materials created # Keys - Material IDs # Values - Materials -OPENMC_MATERIALS = dict() +OPENMC_MATERIALS = {} # A dictionary of all OpenCG Materials created # Keys - Material IDs # Values - Materials -OPENCG_MATERIALS = dict() +OPENCG_MATERIALS = {} # A dictionary of all OpenMC Surfaces created # Keys - Surface IDs # Values - Surfaces -OPENMC_SURFACES = dict() +OPENMC_SURFACES = {} # A dictionary of all OpenCG Surfaces created # Keys - Surface IDs # Values - Surfaces -OPENCG_SURFACES = dict() +OPENCG_SURFACES = {} # A dictionary of all OpenMC Cells created # Keys - Cell IDs # Values - Cells -OPENMC_CELLS = dict() +OPENMC_CELLS = {} # A dictionary of all OpenCG Cells created # Keys - Cell IDs # Values - Cells -OPENCG_CELLS = dict() +OPENCG_CELLS = {} # A dictionary of all OpenMC Universes created # Keys - Universes IDs # Values - Universes -OPENMC_UNIVERSES = dict() +OPENMC_UNIVERSES = {} # A dictionary of all OpenCG Universes created # Keys - Universes IDs # Values - Universes -OPENCG_UNIVERSES = dict() +OPENCG_UNIVERSES = {} # A dictionary of all OpenMC Lattices created # Keys - Lattice IDs # Values - Lattices -OPENMC_LATTICES = dict() +OPENMC_LATTICES = {} # A dictionary of all OpenCG Lattices created # Keys - Lattice IDs # Values - Lattices -OPENCG_LATTICES = dict() +OPENCG_LATTICES = {} @@ -408,7 +408,7 @@ def get_compatible_opencg_cells(opencg_cell, opencg_surface, halfspace): raise ValueError(msg) # Initialize an empty list for the new compatible cells - compatible_cells = list() + compatible_cells = [] # SquarePrism Surfaces if opencg_surface._type in ['x-squareprism', diff --git a/src/utils/openmc/plots.py b/src/utils/openmc/plots.py index dcb74c8adc..dbea64b09e 100644 --- a/src/utils/openmc/plots.py +++ b/src/utils/openmc/plots.py @@ -398,7 +398,7 @@ class PlotsFile(object): def __init__(self): # Initialize PlotsFile class attributes - self._plots = list() + self._plots = [] self._plots_file = ET.Element("plots") diff --git a/src/utils/openmc/statepoint.py b/src/utils/openmc/statepoint.py index dad2f4d268..4ded8d1b65 100644 --- a/src/utils/openmc/statepoint.py +++ b/src/utils/openmc/statepoint.py @@ -186,7 +186,7 @@ class StatePoint(object): # Initialize dictionaries for the Meshes # Keys - Mesh IDs # Values - Mesh objects - self._meshes = dict() + self._meshes = {} # Read the number of Meshes self._n_meshes = self._get_int(path='tallies/meshes/n_meshes')[0] @@ -203,8 +203,8 @@ class StatePoint(object): path='tallies/meshes/keys') else: - self._mesh_keys = list() - self._mesh_ids = list() + self._mesh_keys = [] + self._mesh_ids = [] # Build dictionary of Meshes base = 'tallies/meshes/mesh ' @@ -251,7 +251,7 @@ class StatePoint(object): # Initialize dictionaries for the Tallies # Keys - Tally IDs # Values - Tally objects - self._tallies = dict() + self._tallies = {} # Read the number of tallies self._n_tallies = self._get_int(path='/tallies/n_tallies')[0] @@ -268,8 +268,8 @@ class StatePoint(object): self._n_tallies, path='tallies/keys') else: - self._tally_keys = list() - self._tally_ids = list() + self._tally_keys = [] + self._tally_ids = [] base = 'tallies/tally ' @@ -376,7 +376,7 @@ class StatePoint(object): path='{0}{1}/n_user_score_bins'.format(base, tally_key))[0] # Read scattering moment order strings (e.g., P3, Y-1,2, etc.) - moments = list() + moments = [] subbase = '{0}{1}/moments/'.format(base, tally_key) # Extract the moment order string for each score @@ -664,25 +664,25 @@ class StatePoint(object): for filter in tally._filters: if filter._type == 'surface': - surface_ids = list() + surface_ids = [] for bin in filter._bins: surface_ids.append(summary.surfaces[bin]._id) filter.set_bin_edges(surface_ids) if filter._type in ['cell', 'distribcell']: - distribcell_ids = list() + distribcell_ids = [] for bin in filter._bins: distribcell_ids.append(summary.cells[bin]._id) filter.set_bin_edges(distribcell_ids) if filter._type == 'universe': - universe_ids = list() + universe_ids = [] for bin in filter._bins: universe_ids.append(summary.universes[bin]._id) filter.set_bin_edges(universe_ids) if filter._type == 'material': - material_ids = list() + material_ids = [] for bin in filter._bins: material_ids.append(summary.materials[bin]._id) filter.set_bin_edges(material_ids) diff --git a/src/utils/openmc/summary.py b/src/utils/openmc/summary.py index 14c06f0460..fe179a971e 100644 --- a/src/utils/openmc/summary.py +++ b/src/utils/openmc/summary.py @@ -63,7 +63,7 @@ class Summary(object): # Initialize dictionary for each Nuclide # Keys - Nuclide ZAIDs # Values - Nuclide objects - self.nuclides = dict() + self.nuclides = {} for key in self._f['nuclides'].keys(): @@ -95,7 +95,7 @@ class Summary(object): # Initialize dictionary for each Material # Keys - Material keys # Values - Material objects - self.materials = dict() + self.materials = {} for key in self._f['materials'].keys(): @@ -109,8 +109,8 @@ class Summary(object): nuclides = self._f['materials'][key]['nuclides'][...] n_sab = self._f['materials'][key]['n_sab'][0] - sab_names = list() - sab_xs = list() + sab_names = [] + sab_xs = [] # Read the names of the S(a,b) tables for this Material for i in range(1, n_sab+1): @@ -154,7 +154,7 @@ class Summary(object): # Initialize dictionary for each Surface # Keys - Surface keys # Values - Surfacee objects - self.surfaces = dict() + self.surfaces = {} for key in self._f['geometry/surfaces'].keys(): @@ -237,7 +237,7 @@ class Summary(object): # Initialize dictionary for each Cell # Keys - Cell keys # Values - Cell objects - self.cells = dict() + self.cells = {} # Initialize dictionary for each Cell's fill # (e.g., Material, Universe or Lattice ID) @@ -245,7 +245,7 @@ class Summary(object): # the corresponding objects # Keys - Cell keys # Values - Filling Material, Universe or Lattice ID - self._cell_fills = dict() + self._cell_fills = {} for key in self._f['geometry/cells'].keys(): @@ -266,7 +266,7 @@ class Summary(object): if 'surfaces' in self._f['geometry/cells'][key].keys(): surfaces = self._f['geometry/cells'][key]['surfaces'][...] else: - surfaces = list() + surfaces = [] # Create this Cell cell = openmc.Cell(cell_id=cell_id) @@ -308,7 +308,7 @@ class Summary(object): # Initialize dictionary for each Universe # Keys - Universe keys # Values - Universe objects - self.universes = dict() + self.universes = {} for key in self._f['geometry/universes'].keys(): @@ -338,7 +338,7 @@ class Summary(object): # Initialize lattices for each Lattice # Keys - Lattice keys # Values - Lattice objects - self.lattices = dict() + self.lattices = {} for key in self._f['geometry/lattices'].keys(): diff --git a/src/utils/openmc/surface.py b/src/utils/openmc/surface.py index 087d3981b3..9e4660eced 100644 --- a/src/utils/openmc/surface.py +++ b/src/utils/openmc/surface.py @@ -25,11 +25,11 @@ class Surface(object): # A dictionary of the quadratic surface coefficients # Key - coefficeint name # Value - coefficient value - self._coeffs = dict() + self._coeffs = {} # An ordered list of the coefficient names to export to XML in the # proper order - self._coeff_keys = list() + self._coeff_keys = [] self.set_id(surface_id) self.set_boundary_type(bc_type) diff --git a/src/utils/openmc/tallies.py b/src/utils/openmc/tallies.py index 36880236cf..4250aa4635 100644 --- a/src/utils/openmc/tallies.py +++ b/src/utils/openmc/tallies.py @@ -62,7 +62,7 @@ class Filter(object): def __hash__(self): - hashable = list() + hashable = [] hashable.append(self._type) hashable.append(self._bins) return hash(tuple(hashable)) @@ -530,9 +530,9 @@ class Tally(object): # Initialize Tally class attributes self._id = None self._label = None - self._filters = list() - self._nuclides = list() - self._scores = list() + self._filters = [] + self._nuclides = [] + self._scores = [] self._estimator = None self._num_score_bins = 0 @@ -565,15 +565,15 @@ class Tally(object): clone._mean = copy.deepcopy(self._mean, memo) clone._std_dev = copy.deepcopy(self._std_dev, memo) - clone._filters = list() + clone._filters = [] for filter in self._filters: clone.add_filter(copy.deepcopy(filter, memo)) - clone._nuclides = list() + clone._nuclides = [] for nuclide in self._nuclides: clone.add_nuclide(copy.deepcopy(nuclide, memo)) - clone._scores = list() + clone._scores = [] for score in self._scores: clone.add_score(score) @@ -610,7 +610,7 @@ class Tally(object): def __hash__(self): - hashable = list() + hashable = [] for filter in self._filters: hashable.append((filter._type, tuple(filter._bins))) @@ -1103,7 +1103,7 @@ class Tally(object): tally_group.create_dataset('scores', data=np.array(self._scores)) # Add a string array of the nuclides to the HDF5 group - nuclides = list() + nuclides = [] for nuclide in self._nuclides: nuclides.append(nuclide._name) @@ -1138,10 +1138,10 @@ class Tally(object): if os.path.exists(filename) and append: tally_results = pickle.load(file(filename, 'rb')) else: - tally_results = dict() + tally_results = {} # Create a nested dictionary within the file for this particular Tally - tally_results['Tally-{0}'.format(self._id)] = dict() + tally_results['Tally-{0}'.format(self._id)] = {} tally_group = tally_results['Tally-{0}'.format(self._id)] # Add basic Tally data to the nested dictionary @@ -1151,7 +1151,7 @@ class Tally(object): tally_group['scores'] = np.array(self._scores) # Add a string array of the nuclides to the HDF5 group - nuclides = list() + nuclides = [] for nuclide in self._nuclides: nuclides.append(nuclide._name) @@ -1159,7 +1159,7 @@ class Tally(object): tally_group['nuclides']= np.array(nuclides) # Create a nested dictionary for the Filters - tally_group['filters'] = dict() + tally_group['filters'] = {} filter_group = tally_group['filters'] for filter in self._filters: @@ -1180,8 +1180,8 @@ class TalliesFile(object): def __init__(self): # Initialize TalliesFile class attributes - self._tallies = list() - self._meshes = list() + self._tallies = [] + self._meshes = [] self._tallies_file = ET.Element("tallies") diff --git a/src/utils/openmc/universe.py b/src/utils/openmc/universe.py index 4c5b688b0a..676c94587f 100644 --- a/src/utils/openmc/universe.py +++ b/src/utils/openmc/universe.py @@ -31,7 +31,7 @@ class Cell(object): self._name = None self._fill = None self._type = None - self._surfaces = dict() + self._surfaces = {} self._rotation = None self._translation = None self._offset = None @@ -226,7 +226,7 @@ class Cell(object): def get_all_nuclides(self): - nuclides = dict() + nuclides = {} if self._type != 'void': nuclides.update(self._fill.get_all_nuclides()) @@ -236,7 +236,7 @@ class Cell(object): def get_all_cells(self): - cells = dict() + cells = {} if self._type == 'fill' or self._type == 'lattice': cells.update(self._fill.get_all_cells()) @@ -246,7 +246,7 @@ class Cell(object): def get_all_universes(self): - universes = dict() + universes = {} if self._type == 'fill': universes[self._fill._id] = self._fill @@ -380,7 +380,7 @@ class Universe(object): # Keys - Cell IDs # Values - Cells - self._cells = dict() + self._cells = {} # Keys - Cell IDs # Values - Offsets @@ -483,7 +483,7 @@ class Universe(object): def get_all_nuclides(self): - nuclides = dict() + nuclides = {} # Append all Nuclides in each Cell in the Universe to the dictionary for cell_id, cell in self._cells.items(): @@ -494,7 +494,7 @@ class Universe(object): def get_all_cells(self): - cells = dict() + cells = {} # Add this Universe's cells to the dictionary cells.update(self._cells) @@ -511,7 +511,7 @@ class Universe(object): # Get all Cells in this Universe cells = self.get_all_cells() - universes = dict() + universes = {} # Append all Universes containing each Cell to the dictionary for cell_id, cell in cells.items(): @@ -633,7 +633,7 @@ class Lattice(object): def get_unique_universes(self): unique_universes = np.unique(self._universes.ravel()) - universes = dict() + universes = {} for universe in unique_universes: universes[universe._id] = universe @@ -643,7 +643,7 @@ class Lattice(object): def get_all_nuclides(self): - nuclides = dict() + nuclides = {} # Get all unique Universes contained in each of the lattice cells unique_universes = self.get_unique_universes() @@ -657,7 +657,7 @@ class Lattice(object): def get_all_cells(self): - cells = dict() + cells = {} unique_universes = self.get_unique_universes() for universe_id, universe in unique_universes.items(): @@ -670,7 +670,7 @@ class Lattice(object): # Initialize a dictionary of all Universes contained by the Lattice # in each nested Universe level - all_universes = dict() + all_universes = {} # Get all unique Universes contained in each of the lattice cells unique_universes = self.get_unique_universes() From 1c3ac71babfbef768c6eb60dd4d3ff8d77ce035e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Apr 2015 12:35:37 -0500 Subject: [PATCH 23/29] Use triple double quotes instead of triple single quotes --- src/utils/openmc/clean_xml.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utils/openmc/clean_xml.py b/src/utils/openmc/clean_xml.py index 437c2df15c..619475fa07 100644 --- a/src/utils/openmc/clean_xml.py +++ b/src/utils/openmc/clean_xml.py @@ -66,11 +66,11 @@ def sort_xml_elements(tree): def clean_xml_indentation(element, level=0): - ''' + """ copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint it basically walks your tree and adds spaces and newlines so the tree is printed in a nice way - ''' + """ i = "\n" + level*" " From db61158325580aaea3c895cd96d5bc04619459ff Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Apr 2015 12:52:00 -0500 Subject: [PATCH 24/29] Make sure /dev/null isn't open for the entire process --- src/utils/openmc/executor.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/utils/openmc/executor.py b/src/utils/openmc/executor.py index f2d26de9b6..abc2e9e49e 100644 --- a/src/utils/openmc/executor.py +++ b/src/utils/openmc/executor.py @@ -3,9 +3,6 @@ import subprocess import os -FNULL = open(os.devnull, 'w') - - class Executor(object): @@ -34,7 +31,8 @@ class Executor(object): subprocess.check_call('openmc -p', shell=True, cwd=self._working_directory) else: - subprocess.check_call('openmc -p', shell=True, stdout=FNULL, + subprocess.check_call('openmc -p', shell=True, + stdout=open(os.devnull, 'w'), cwd=self._working_directory) @@ -69,5 +67,6 @@ class Executor(object): subprocess.check_call(command, shell=True, cwd=self._working_directory) else: - subprocess.check_call(command, shell=True, stdout=FNULL, + subprocess.check_call(command, shell=True, + stdout=open(os.devnull, 'w'), cwd=self._working_directory) From 37a9b7478f1fbe19f55150c96d925a98ab2d871a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Apr 2015 13:13:54 -0500 Subject: [PATCH 25/29] Respond to various scopatz comments in pull request #331 --- src/utils/openmc/cmfd.py | 6 ++- src/utils/openmc/executor.py | 3 +- src/utils/openmc/geometry.py | 4 +- src/utils/openmc/material.py | 16 ++++---- src/utils/openmc/nuclide.py | 5 +-- src/utils/openmc/opencg_compatible.py | 6 ++- src/utils/openmc/plots.py | 6 ++- src/utils/openmc/settings.py | 57 ++++++++++++++------------- src/utils/openmc/statepoint.py | 6 ++- src/utils/openmc/summary.py | 3 +- src/utils/openmc/surface.py | 3 +- src/utils/openmc/tallies.py | 8 ++-- src/utils/openmc/universe.py | 10 +++-- 13 files changed, 74 insertions(+), 59 deletions(-) diff --git a/src/utils/openmc/cmfd.py b/src/utils/openmc/cmfd.py index 443867284a..5d10ba5383 100644 --- a/src/utils/openmc/cmfd.py +++ b/src/utils/openmc/cmfd.py @@ -1,7 +1,9 @@ +from xml.etree import ElementTree as ET + +import numpy as np + from openmc.checkvalue import * from openmc.clean_xml import * -from xml.etree import ElementTree as ET -import numpy as np class CMFDMesh(object): diff --git a/src/utils/openmc/executor.py b/src/utils/openmc/executor.py index abc2e9e49e..446111c3d9 100644 --- a/src/utils/openmc/executor.py +++ b/src/utils/openmc/executor.py @@ -1,7 +1,8 @@ -from openmc.checkvalue import * import subprocess import os +from openmc.checkvalue import * + class Executor(object): diff --git a/src/utils/openmc/geometry.py b/src/utils/openmc/geometry.py index 17d0f5fa87..db2227ba5a 100644 --- a/src/utils/openmc/geometry.py +++ b/src/utils/openmc/geometry.py @@ -1,6 +1,8 @@ +from xml.etree import ElementTree as ET + import openmc from openmc.clean_xml import * -from xml.etree import ElementTree as ET + def reset_auto_ids(): openmc.reset_auto_material_id() diff --git a/src/utils/openmc/material.py b/src/utils/openmc/material.py index 0d6e39b992..99402ef75e 100644 --- a/src/utils/openmc/material.py +++ b/src/utils/openmc/material.py @@ -1,12 +1,13 @@ +from collections import MappingView +from copy import deepcopy import warnings +from xml.etree import ElementTree as ET + +import numpy as np import openmc from openmc.checkvalue import * from openmc.clean_xml import * -from xml.etree import ElementTree as ET -from collections import MappingView -from copy import deepcopy -import numpy as np # A list of all IDs for all Materials created @@ -70,14 +71,13 @@ class Material(object): def set_id(self, material_id=None): - global MATERIAL_IDS + global AUTO_MATERIAL_ID, MATERIAL_IDS # If the Material already has an ID, remove it from global list if not self._id is None: MATERIAL_IDS.remove(self._id) if material_id is None: - global AUTO_MATERIAL_ID self._id = AUTO_MATERIAL_ID MATERIAL_IDS.append(AUTO_MATERIAL_ID) AUTO_MATERIAL_ID += 1 @@ -260,13 +260,13 @@ class Material(object): string += '{0: <16}{1}{2}'.format('\tDensity', '=\t', self._density) string += ' [{0}]\n'.format(self._density_units) - string += '{0: <16}'.format('\tS(a,b) Tables') + '\n' + string += '{0: <16}\n'.format('\tS(a,b) Tables') for sab in self._sab: string += '{0: <16}{1}[{2}{3}]\n'.format('\tS(a,b)', '=\t', sab[0], sab[1]) - string += '{0: <16}'.format('\tNuclides') + '\n' + string += '{0: <16}\n'.format('\tNuclides') for nuclide in self._nuclides: percent = self._nuclides[nuclide][1] diff --git a/src/utils/openmc/nuclide.py b/src/utils/openmc/nuclide.py index e100d6d74e..5281eacd4c 100644 --- a/src/utils/openmc/nuclide.py +++ b/src/utils/openmc/nuclide.py @@ -36,10 +36,7 @@ class Nuclide(object): def __hash__(self): - hashable = [] - hashable.append(self._name) - hashable.append(self._xs) - return hash(tuple(hashable)) + return hash((self._name, self._xs)) def set_name(self, name): diff --git a/src/utils/openmc/opencg_compatible.py b/src/utils/openmc/opencg_compatible.py index f80eb2c2c7..b7644aa48c 100644 --- a/src/utils/openmc/opencg_compatible.py +++ b/src/utils/openmc/opencg_compatible.py @@ -1,7 +1,9 @@ -import openmc -import opencg import copy + import numpy as np +import opencg + +import openmc # A dictionary of all OpenMC Materials created diff --git a/src/utils/openmc/plots.py b/src/utils/openmc/plots.py index dbea64b09e..a69fb1abc1 100644 --- a/src/utils/openmc/plots.py +++ b/src/utils/openmc/plots.py @@ -1,7 +1,9 @@ +from xml.etree import ElementTree as ET + +import numpy as np + from openmc.checkvalue import * from openmc.clean_xml import * -from xml.etree import ElementTree as ET -import numpy as np # A static variable for auto-generated Plot IDs diff --git a/src/utils/openmc/settings.py b/src/utils/openmc/settings.py index 0e314545a2..240892e420 100644 --- a/src/utils/openmc/settings.py +++ b/src/utils/openmc/settings.py @@ -1,10 +1,11 @@ +import collections import warnings +from xml.etree import ElementTree as ET + +import numpy as np from openmc.checkvalue import * from openmc.clean_xml import * -from xml.etree import ElementTree as ET -import numpy as np -import multiprocessing class SettingsFile(object): @@ -152,16 +153,16 @@ class SettingsFile(object): self._source_file = source_file - def set_source_space(self, type, params): + def set_source_space(self, stype, params): - if not is_string(type): + if not is_string(stype): msg = 'Unable to set source space type to a non-string ' \ - 'value {0}'.format(type) + 'value {0}'.format(stype) raise ValueError(msg) - elif not type in ['box', 'point']: + elif not stype in ['box', 'point']: msg = 'Unable to set source space type to {0} since it is not ' \ - 'box or point'.format(type) + 'box or point'.format(stype) raise ValueError(msg) elif not isinstance(params, (tuple, list, np.ndarray)): @@ -181,20 +182,20 @@ class SettingsFile(object): 'is not an integer or floating point value'.format(param) raise ValueError(msg) - self._source_space_type = type + self._source_space_type = stype self._source_space_params = params - def set_source_angle(self, type, params=[]): + def set_source_angle(self, stype, params=[]): - if not is_string(type): + if not is_string(stype): msg = 'Unable to set source angle type to a non-string ' \ - 'value {0}'.format(type) + 'value {0}'.format(stype) raise ValueError(msg) - elif not type in ['isotropic', 'monodirectional']: + elif not stype in ['isotropic', 'monodirectional']: msg = 'Unable to set source angle type to {0} since it is not ' \ - 'isotropic or monodirectional'.format(type) + 'isotropic or monodirectional'.format(stype) raise ValueError(msg) elif not isinstance(params, (tuple, list, np.ndarray)): @@ -202,12 +203,12 @@ class SettingsFile(object): 'not a Python list/tuple or NumPy array'.format(params) raise ValueError(msg) - elif type is 'isotropic' and not params is None: + elif stype == 'isotropic' and not params is None: msg = 'Unable to set source angle parameters since they are not ' \ 'it is not supported for isotropic type sources' raise ValueError(msg) - elif type is 'monodirectional' and len(params) != 3: + elif stype == 'monodirectional' and len(params) != 3: msg = 'Unable to set source angle parameters to {0} ' \ 'since 3 parameters are required for monodirectional ' \ 'sources'.format(params) @@ -220,20 +221,20 @@ class SettingsFile(object): 'is not an integer or floating point value'.format(param) raise ValueError(msg) - self._source_angle_type = type + self._source_angle_type = stype self._source_angle_params = params - def set_source_energy(self, type, params=[]): + def set_source_energy(self, stype, params=[]): - if not is_string(type): + if not is_string(stype): msg = 'Unable to set source energy type to a non-string ' \ - 'value {0}'.format(type) + 'value {0}'.format(stype) raise ValueError(msg) - elif not type in ['monoenergetic', 'watt', 'maxwell']: + elif not stype in ['monoenergetic', 'watt', 'maxwell']: msg = 'Unable to set source energy type to {0} since it is not ' \ - 'monoenergetic, watt or maxwell'.format(type) + 'monoenergetic, watt or maxwell'.format(stype) raise ValueError(msg) elif not isinstance(params, (tuple, list, np.ndarray)): @@ -241,19 +242,19 @@ class SettingsFile(object): 'is not a Python list/tuple or NumPy array'.format(params) raise ValueError(msg) - elif type is 'monoenergetic' and not len(params) != 1: + elif stype == 'monoenergetic' and not len(params) != 1: msg = 'Unable to set source energy parameters to {0} ' \ 'since 1 paramater is required for monenergetic ' \ 'sources'.format(params) raise ValueError(msg) - elif type is 'watt' and len(params) != 2: + elif stype == 'watt' and len(params) != 2: msg = 'Unable to set source energy parameters to {0} ' \ 'since 2 parameters are required for monoenergetic ' \ 'sources'.format(params) raise ValueError(msg) - elif type is 'maxwell' and len(params) != 2: + elif stype == 'maxwell' and len(params) != 2: msg = 'Unable to set source energy parameters to {0} since 1 ' \ 'parameter is required for maxwell sources'.format(params) raise ValueError(msg) @@ -266,7 +267,7 @@ class SettingsFile(object): 'value'.format(param) raise ValueError(msg) - self._source_energy_type = type + self._source_energy_type = stype self._source_energy_params = params @@ -806,7 +807,7 @@ class SettingsFile(object): warnings.warn('This feature is not yet implemented in a release ' \ 'version of openmc') - if not type(allow) == bool: + if not isinstance(allow, bool): msg = 'Unable to set DD allow_leakage {0} which is ' \ 'not a Python bool'.format(dimension) raise ValueError(msg) @@ -820,7 +821,7 @@ class SettingsFile(object): warnings.warn('This feature is not yet implemented in a release ' \ 'version of openmc') - if not type(interactions) == bool: + if not isinstance(interactions, bool): msg = 'Unable to set DD count_interactions {0} which is ' \ 'not a Python bool'.format(dimension) raise ValueError(msg) diff --git a/src/utils/openmc/statepoint.py b/src/utils/openmc/statepoint.py index 4ded8d1b65..eadfcc4c1f 100644 --- a/src/utils/openmc/statepoint.py +++ b/src/utils/openmc/statepoint.py @@ -1,11 +1,13 @@ -import struct, copy +import copy +import struct + import numpy as np import scipy.stats + import openmc from openmc.constants import * - class SourceSite(object): def __init__(self): diff --git a/src/utils/openmc/summary.py b/src/utils/openmc/summary.py index fe179a971e..caddc99a98 100644 --- a/src/utils/openmc/summary.py +++ b/src/utils/openmc/summary.py @@ -1,6 +1,7 @@ +import numpy as np + import openmc from openmc.opencg_compatible import get_opencg_geometry -import numpy as np try: import h5py diff --git a/src/utils/openmc/surface.py b/src/utils/openmc/surface.py index 9e4660eced..2ed28e3835 100644 --- a/src/utils/openmc/surface.py +++ b/src/utils/openmc/surface.py @@ -1,6 +1,7 @@ +from xml.etree import ElementTree as ET + from openmc.checkvalue import * from openmc.constants import BC_TYPES -from xml.etree import ElementTree as ET # A static variable for auto-generated Surface IDs diff --git a/src/utils/openmc/tallies.py b/src/utils/openmc/tallies.py index 4250aa4635..9d9f8acaf2 100644 --- a/src/utils/openmc/tallies.py +++ b/src/utils/openmc/tallies.py @@ -1,10 +1,12 @@ +from xml.etree import ElementTree as ET +import os, copy + +import numpy as np + from openmc import Nuclide from openmc.clean_xml import * from openmc.checkvalue import * from openmc.constants import * -from xml.etree import ElementTree as ET -import numpy as np -import os, copy # "Static" variables for auto-generated Tally and Mesh IDs diff --git a/src/utils/openmc/universe.py b/src/utils/openmc/universe.py index 676c94587f..5bb9c29304 100644 --- a/src/utils/openmc/universe.py +++ b/src/utils/openmc/universe.py @@ -1,9 +1,11 @@ +import abc +from collections import OrderedDict +from xml.etree import ElementTree as ET + +import numpy as np + import openmc from openmc.checkvalue import * -from xml.etree import ElementTree as ET -from collections import OrderedDict -import numpy as np -import abc ################################################################################ From 4bf5f3918ee0d18d202c550e52e470bcafb3127f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Apr 2015 14:22:03 -0500 Subject: [PATCH 26/29] Remove two imports on one line --- src/utils/openmc/tallies.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/utils/openmc/tallies.py b/src/utils/openmc/tallies.py index 9d9f8acaf2..1f7271de81 100644 --- a/src/utils/openmc/tallies.py +++ b/src/utils/openmc/tallies.py @@ -1,5 +1,6 @@ +import copy +import os from xml.etree import ElementTree as ET -import os, copy import numpy as np From 52a4a10d54374b08c5835d02077fd1edcdc547ac Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Apr 2015 14:35:49 -0500 Subject: [PATCH 27/29] Fix import for statepoint for test_union_energy_grids --- tests/test_union_energy_grids/results.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/test_union_energy_grids/results.py b/tests/test_union_energy_grids/results.py index be13ee66f1..b616e2874b 100644 --- a/tests/test_union_energy_grids/results.py +++ b/tests/test_union_energy_grids/results.py @@ -2,20 +2,20 @@ import sys -# import statepoint sys.path.insert(0, '../../src/utils') -import statepoint +from openmc.statepoint import StatePoint # read in statepoint file if len(sys.argv) > 1: - sp = statepoint.StatePoint(sys.argv[1]) + print(sys.argv) + sp = StatePoint(sys.argv[1]) else: - sp = statepoint.StatePoint('statepoint.10.binary') + sp = 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]) From 63f8535e840adb02a62e3025ada46f60c5e990a1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Apr 2015 14:39:43 -0500 Subject: [PATCH 28/29] Get rid of unused variables in state_point.F90. Also remove unnecessary null() initializations. --- src/state_point.F90 | 34 +++++++++++++--------------------- 1 file changed, 13 insertions(+), 21 deletions(-) diff --git a/src/state_point.F90 b/src/state_point.F90 index c3edc3bff9..e30be0ea31 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -39,17 +39,13 @@ contains subroutine write_state_point() character(MAX_FILE_LEN) :: filename - integer :: i, j, k, m - integer :: n_x, n_y, n_z - character(11), allocatable :: name_array(:) + integer :: i, j, k integer, allocatable :: id_array(:) integer, allocatable :: key_array(:) - type(StructuredMesh), pointer :: mesh => null() - type(TallyObject), pointer :: tally => null() - type(ElemKeyValueII), pointer :: current => null() - type(ElemKeyValueII), pointer :: next => null() - type(ElemKeyValueCI), pointer :: cur_nuclide => null() - type(ElemKeyValueCI), pointer :: next_nuclide => null() + type(StructuredMesh), pointer :: mesh + type(TallyObject), pointer :: tally + type(ElemKeyValueII), pointer :: current + type(ElemKeyValueII), pointer :: next character(8) :: moment_name ! name of moment (e.g, P3) integer :: n_order ! loop index for moment orders integer :: nm_order ! loop index for Ynm moment orders @@ -306,14 +302,14 @@ contains case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) moment_name = 'P' // to_str(tally % moment_order(k)) call sp % write_data(moment_name, "order" // trim(to_str(k)), & - group="tallies/tally " // trim(to_str(tally % id)) // & + group="tallies/tally " // trim(to_str(tally % id)) // & "/moments") k = k + 1 case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN) do n_order = 0, tally % moment_order(k) moment_name = 'P' // trim(to_str(n_order)) call sp % write_data(moment_name, "order" // trim(to_str(k)), & - group="tallies/tally " // trim(to_str(tally % id)) // & + group="tallies/tally " // trim(to_str(tally % id)) // & "/moments") k = k + 1 end do @@ -325,7 +321,7 @@ contains trim(to_str(nm_order)) call sp % write_data(moment_name, "order" // & trim(to_str(k)), & - group="tallies/tally " // trim(to_str(tally % id)) // & + group="tallies/tally " // trim(to_str(tally % id)) // & "/moments") k = k + 1 end do @@ -504,9 +500,9 @@ contains real(8) :: dummy ! temporary receive buffer for non-root reduces #endif integer, allocatable :: id_array(:) - type(ElemKeyValueII), pointer :: current => null() - type(ElemKeyValueII), pointer :: next => null() - type(TallyObject), pointer :: tally => null() + type(ElemKeyValueII), pointer :: current + type(ElemKeyValueII), pointer :: next + type(TallyObject), pointer :: tally type(TallyResult), allocatable :: tallyresult_temp(:,:) ! ========================================================================== @@ -661,14 +657,10 @@ contains integer, allocatable :: key_array(:) integer :: curr_key integer, allocatable :: temp_array(:) - integer, allocatable :: temp_array3D(:,:,:) - integer, allocatable :: temp_array4D(:,:,:,:) logical :: source_present - real(8) :: l real(8) :: real_array(3) - real(8), allocatable :: temp_real_array(:) - type(StructuredMesh), pointer :: mesh => null() - type(TallyObject), pointer :: tally => null() + type(StructuredMesh), pointer :: mesh + type(TallyObject), pointer :: tally integer :: n_order ! loop index for moment orders integer :: nm_order ! loop index for Ynm moment orders character(8) :: moment_name ! name of moment (e.g, P3, Y-1,1) From ee5271e0e3f974f66749601f5de488db5541a0de Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Apr 2015 14:44:40 -0500 Subject: [PATCH 29/29] Remove unused hdf5_read_string_1darray routine --- src/hdf5_interface.F90 | 71 +++++++++++------------------------------- 1 file changed, 18 insertions(+), 53 deletions(-) diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index aedd53dd19..2252c3e763 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -68,7 +68,6 @@ module hdf5_interface module procedure hdf5_read_integer_4Darray module procedure hdf5_read_long module procedure hdf5_read_string - module procedure hdf5_read_string_1Darray #ifdef MPI module procedure hdf5_read_double_parallel module procedure hdf5_read_double_1Darray_parallel @@ -82,7 +81,6 @@ module hdf5_interface module procedure hdf5_read_integer_4Darray_parallel module procedure hdf5_read_long_parallel module procedure hdf5_read_string_parallel -! module procedure hdf5_read_string_1Darray_parallel #endif end interface hdf5_read_data @@ -208,7 +206,7 @@ contains logical :: status ! does the group exist ! Check if group exists - call h5ltpath_valid_f(hdf5_fh, trim(group), .true., status, hdf5_err) + call h5ltpath_valid_f(hdf5_fh, trim(group), .true., status, hdf5_err) ! Either create or open group if (status) then @@ -259,7 +257,7 @@ contains integer(HID_T), intent(in) :: group ! name of group character(*), intent(in) :: name ! name of data - integer, intent(inout) :: buffer ! read data to here + integer, intent(inout) :: buffer ! read data to here integer :: buffer_copy(1) ! need an array for read @@ -463,7 +461,7 @@ contains integer(HID_T), intent(in) :: group ! name of group character(*), intent(in) :: name ! name of data - real(8), intent(inout) :: buffer ! read data to here + real(8), intent(inout) :: buffer ! read data to here real(8) :: buffer_copy(1) ! need an array for read @@ -728,7 +726,7 @@ contains call h5dcreate_f(group, name, H5T_STRING, dspace, dset, hdf5_err) ! Set up dimesnions of string to write - dims2 = (/length, 1/) ! full array of strings to write + dims2 = (/length, 1/) ! full array of strings to write dims1(1) = length ! length of string ! Copy over string buffer to a rank 1 array @@ -762,7 +760,7 @@ contains ! character(len=length, kind=c_char), pointer :: chr_ptr ! f_ptr = c_loc(buf_ptr(1)) ! call h5dread_f(dset, H5T_STRING, f_ptr, hdf5_err, xfer_prp=plist) -! call c_f_pointer(buf_ptr(1), chr_ptr) +! call c_f_pointer(buf_ptr(1), chr_ptr) ! buffer = chr_ptr ! nullify(chr_ptr) @@ -788,39 +786,6 @@ contains end subroutine hdf5_read_string -!=============================================================================== -! HDF5_READ_STRING_1DARRAY reads string data -!=============================================================================== - - subroutine hdf5_read_string_1Darray(group, name, buffer, length, rank) - - integer(HID_T), intent(in) :: group ! name of group - character(*), intent(in) :: name ! name of data - integer, intent(in) :: length ! length of strings - integer, intent(in) :: rank ! number of strings - character(*), intent(inout) :: buffer(length,rank) ! read data to here - - ! Open dataset - call h5dopen_f(group, name, dset, hdf5_err) - - ! Get dataspace to read - call h5dget_space_f(dset, dspace, hdf5_err) - - ! Set dimensions - dims2 = (/length, rank/) - dims1(1) = length - - ! Read in the data - !call h5ltread_dataset_string_f(dset, H5T_STRING, buffer, dims2, dims1, hdf5_err, & - ! mem_space_id=dspace, xfer_prp = plist) - - !call h5ltread_dataset_f(group, name, dims, dims2, H5T_STRING, buffer, hdf5_err) - - ! Close dataset - call h5dclose_f(dset, hdf5_err) - - end subroutine hdf5_read_string_1Darray - !=============================================================================== ! HDF5_WRITE_ATTRIBUTE_STRING writes a string attribute to a variables !=============================================================================== @@ -873,7 +838,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -951,7 +916,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1031,7 +996,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1112,7 +1077,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1194,7 +1159,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1273,7 +1238,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1351,7 +1316,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1431,7 +1396,7 @@ contains f_ptr = c_loc(buffer(1,1)) call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1512,7 +1477,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1594,7 +1559,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1674,7 +1639,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, long_type, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1761,7 +1726,7 @@ contains call h5dcreate_f(group, name, H5T_STRING, dspace, dset, hdf5_err) ! Set up dimesnions of string to write - dims2 = (/length, 1/) ! full array of strings to write + dims2 = (/length, 1/) ! full array of strings to write dims1(1) = length ! length of string ! Copy over string buffer to a rank 1 array @@ -1797,7 +1762,7 @@ contains ! character(len=length, kind=c_char), pointer :: chr_ptr ! f_ptr = c_loc(buf_ptr(1)) ! call h5dread_f(dset, H5T_STRING, f_ptr, hdf5_err, xfer_prp=plist) -! call c_f_pointer(buf_ptr(1), chr_ptr) +! call c_f_pointer(buf_ptr(1), chr_ptr) ! buffer = chr_ptr ! nullify(chr_ptr)