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

Merge pull request #367 from walshjon/mat-union

Browse source 
unionized material energy-cross section grid lookup
This commit is contained in:
Paul Romano 2015-04-16 08:13:49 -05:00
commit 4677f2fb40
17 changed files with 376 additions and 35 deletions
Download patch file Download diff file Expand all files Collapse all files

14
docs/source/usersguide/input.rst
View file

@ -174,11 +174,15 @@ should be performed. It has the following attributes/sub-elements:
-------------------------
The ``<energy_grid>`` 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

5
src/ace_header.F90
View file

@ -103,7 +103,7 @@ module ace_header
! Energy grid information
integer :: n_grid ! # of nuclide grid points
integer, allocatable :: grid_index(:) ! union grid pointers / log grid mapping
integer, allocatable :: grid_index(:) ! log grid mapping indices
real(8), allocatable :: energy(:) ! energy values corresponding to xs
! Microscopic cross sections
@ -372,9 +372,6 @@ module ace_header
integer :: i ! Loop counter
if (allocated(this % grid_index)) &
deallocate(this % grid_index)
if (allocated(this % energy)) &
deallocate(this % energy, this % total, this % elastic, &
& this % fission, this % nu_fission, this % absorption)

5
src/constants.F90
View file

@ -368,8 +368,9 @@ module constants
! Energy grid methods
integer, parameter :: &
GRID_NUCLIDE = 1, & ! non-unionized energy grid
GRID_LOGARITHM = 2 ! logarithmic mapping
GRID_NUCLIDE = 1, & ! unique energy grid for each nuclide
GRID_MAT_UNION = 2, & ! material union grids with pointers
GRID_LOGARITHM = 3 ! lethargy mapping
! Running modes
integer, parameter :: &

75
src/cross_section.F90
View file

@ -15,6 +15,9 @@ module cross_section
implicit none
save
integer :: union_grid_index
!$omp threadprivate(union_grid_index)
contains
!===============================================================================
@ -36,11 +39,11 @@ contains
!$omp threadprivate(mat)
! Set all material macroscopic cross sections to zero
material_xs % total = ZERO
material_xs % elastic = ZERO
material_xs % absorption = ZERO
material_xs % fission = ZERO
material_xs % nu_fission = ZERO
material_xs % total = ZERO
material_xs % elastic = ZERO
material_xs % absorption = ZERO
material_xs % fission = ZERO
material_xs % nu_fission = ZERO
material_xs % kappa_fission = ZERO
! Exit subroutine if material is void
@ -48,6 +51,10 @@ contains
mat => materials(p % material)
! Find energy index on global or material unionized grid
if (grid_method == GRID_MAT_UNION) &
call find_energy_index(p % E, p % material)
! Determine if this material has S(a,b) tables
check_sab = (mat % n_sab > 0)
@ -88,9 +95,9 @@ contains
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E)
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i)
else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E)
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i)
end if
! ========================================================================
@ -131,24 +138,32 @@ contains
! given index in the nuclides array at the energy of the given particle
!===============================================================================
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E)
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
integer :: i_grid ! index on nuclide energy grid
integer :: i_low, i_high ! bounding indices from logarithmic mapping
integer :: u ! index into logarithmic mapping array
integer, intent(in) :: i_mat ! index into materials array
integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material
integer :: i_grid ! index on nuclide energy grid
integer :: i_low ! lower logarithmic mapping index
integer :: i_high ! upper logarithmic mapping index
integer :: u ! index into logarithmic mapping array
real(8), intent(in) :: E ! energy
real(8) :: f ! interp factor on nuclide energy grid
type(Nuclide), pointer, save :: nuc => null()
!$omp threadprivate(nuc)
type(Nuclide), pointer, save :: nuc => null()
type(Material), pointer, save :: mat => null()
!$omp threadprivate(nuc, mat)
! Set pointer to nuclide
! Set pointer to nuclide and material
nuc => nuclides(i_nuclide)
mat => materials(i_mat)
! Determine index on nuclide energy grid
select case (grid_method)
case (GRID_MAT_UNION)
i_grid = mat % nuclide_grid_index(i_nuc_mat, union_grid_index)
case (GRID_LOGARITHM)
! Determine the energy grid index using a logarithmic mapping to reduce
! the energy range over which a binary search needs to be performed
@ -173,7 +188,7 @@ contains
! Perform binary search on the nuclide energy grid in order to determine
! which points to interpolate between
if (E < nuc % energy(1)) then
if (E <= nuc % energy(1)) then
i_grid = 1
elseif (E > nuc % energy(nuc % n_grid)) then
i_grid = nuc % n_grid - 1
@ -506,6 +521,32 @@ contains
end subroutine calculate_urr_xs
!===============================================================================
! FIND_ENERGY_INDEX determines the index on the union energy grid at a certain
! energy
!===============================================================================
subroutine find_energy_index(E, i_mat)
real(8), intent(in) :: E ! energy of particle
integer, intent(in) :: i_mat ! material index
type(Material), pointer, save :: mat => null() ! pointer to current material
!$omp threadprivate(mat)
mat => materials(i_mat)
! if the energy is outside of energy grid range, set to first or last
! index. Otherwise, do a binary search through the union energy grid.
if (E <= mat % e_grid(1)) then
union_grid_index = 1
elseif (E > mat % e_grid(mat % n_grid)) then
union_grid_index = mat % n_grid - 1
else
union_grid_index = binary_search(mat % e_grid, mat % n_grid, E)
end if
end subroutine find_energy_index
!===============================================================================
! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section
! for a given nuclide at the trial relative energy used in resonance scattering

155
src/energy_grid.F90
View file

@ -1,6 +1,9 @@
module energy_grid
use global
use list_header, only: ListReal
use output, only: write_message
implicit none
@ -10,6 +13,53 @@ module energy_grid
contains
!===============================================================================
! UNIONIZED_GRID creates a unionized energy grid, for the entire problem or for
! each material, composed of the grids from each nuclide in the entire problem,
! or each material, respectively. Right now, the grid for each nuclide is added
! into a linked list one at a time with an effective insertion sort. Could be
! done with a hash for all energy points and then a quicksort at the end (what
! hash function to use?)
!===============================================================================
subroutine unionized_grid()
integer :: i ! index in nuclides array
integer :: j ! index in materials array
type(ListReal), pointer, save :: list => null()
type(Nuclide), pointer, save :: nuc => null()
type(Material), pointer, save :: mat => null()
!$omp threadprivate(list, nuc, mat)
call write_message("Creating unionized energy grid...", 5)
! add grid points for each nuclide in the material
do j = 1, n_materials
mat => materials(j)
do i = 1, mat % n_nuclides
nuc => nuclides(mat % nuclide(i))
call add_grid_points(list, nuc % energy)
end do
! set size of unionized material energy grid
mat % n_grid = list % size()
! create allocated array from linked list
allocate(mat % e_grid(mat % n_grid))
do i = 1, mat % n_grid
mat % e_grid(i) = list % get_item(i)
end do
! delete linked list and dictionary
call list % clear()
deallocate(list)
end do
! Set pointers to unionized energy grid for each nuclide
call grid_pointers()
end subroutine unionized_grid
!===============================================================================
! LOGARITHMIC_GRID determines a logarithmic mapping for energies to bounding
! indices on a nuclide energy grid
@ -59,4 +109,109 @@ contains
end subroutine logarithmic_grid
!===============================================================================
! ADD_GRID_POINTS adds energy points from the 'energy' array into a linked list
! of points already stored from previous arrays.
!===============================================================================
subroutine add_grid_points(list, energy)
type(ListReal), pointer :: list
real(8), intent(in) :: energy(:)
integer :: i ! index in energy array
integer :: n ! size of energy array
integer :: current ! current index
real(8) :: E ! actual energy value
i = 1
n = size(energy)
! If the original list is empty, we need to allocate the first element and
! store first energy point
if (.not. associated(list)) then
allocate(list)
do i = 1, n
call list % append(energy(i))
end do
return
end if
! Set current index to beginning of the list
current = 1
do while (i <= n)
E = energy(i)
! If we've reached the end of the grid energy list, add the remaining
! energy points to the end
if (current > list % size()) then
! Finish remaining energies
do while (i <= n)
call list % append(energy(i))
i = i + 1
end do
exit
end if
if (E < list % get_item(current)) then
! Insert new energy in this position
call list % insert(current, E)
! Advance index in linked list and in new energy grid
i = i + 1
current = current + 1
elseif (E == list % get_item(current)) then
! Found the exact same energy, no need to store duplicates so just
! skip and move to next index
i = i + 1
current = current + 1
else
current = current + 1
end if
end do
end subroutine add_grid_points
!===============================================================================
! GRID_POINTERS creates an array of pointers (ints) for each nuclide to link
! each point on the nuclide energy grid to one on a unionized energy grid
!===============================================================================
subroutine grid_pointers()
integer :: i ! loop index for nuclides
integer :: j ! loop index for nuclide energy grid
integer :: k ! loop index for materials
integer :: index_e ! index on union energy grid
real(8) :: union_energy ! energy on union grid
real(8) :: energy ! energy on nuclide grid
type(Nuclide), pointer :: nuc => null()
type(Material), pointer :: mat => null()
do k = 1, n_materials
mat => materials(k)
allocate(mat % nuclide_grid_index(mat % n_nuclides, mat % n_grid))
do i = 1, mat % n_nuclides
nuc => nuclides(mat % nuclide(i))
index_e = 1
energy = nuc % energy(index_e)
do j = 1, mat % n_grid
union_energy = mat % e_grid(j)
if (union_energy >= energy .and. index_e < nuc % n_grid) then
index_e = index_e + 1
energy = nuc % energy(index_e)
end if
mat % nuclide_grid_index(i,j) = index_e - 1
end do
end do
end do
end subroutine grid_pointers
end module energy_grid

1
src/global.F90
View file

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

15
src/initialize.F90
View file

@ -4,7 +4,7 @@ module initialize
use bank_header, only: Bank
use constants
use dict_header, only: DictIntInt, ElemKeyValueII
use energy_grid, only: logarithmic_grid, grid_method
use energy_grid, only: logarithmic_grid, grid_method, unionized_grid
use error, only: fatal_error, warning
use geometry, only: neighbor_lists
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
@ -109,10 +109,17 @@ contains
! Create linked lists for multiple instances of the same nuclide
call same_nuclide_list()
! Construct logarithmic energy grid for cross-sections
if (grid_method == GRID_LOGARITHM) then
! Construct unionized or log energy grid for cross-sections
select case (grid_method)
case (GRID_NUCLIDE)
continue
case (GRID_MAT_UNION)
call time_unionize % start()
call unionized_grid()
call time_unionize % stop()
case (GRID_LOGARITHM)
call logarithmic_grid()
end if
end select
! Allocate and setup tally stride, matching_bins, and tally maps
call configure_tallies()

7
src/input_xml.F90
View file

@ -213,8 +213,11 @@ contains
select case (trim(temp_str))
case ('nuclide')
grid_method = GRID_NUCLIDE
case ('union')
call fatal_error("Union energy grid is no longer supported.")
case ('material-union', 'union')
grid_method = GRID_MAT_UNION
if (trim(temp_str) == 'union') &
call warning('Energy grids will be unionized by material. Global&
& energy grid unionization is no longer an allowed option.')
case ('logarithm', 'logarithmic', 'log')
grid_method = GRID_LOGARITHM
case default

7
src/material_header.F90
View file

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

2
src/relaxng/settings.rnc
View file

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

2
src/relaxng/settings.rng
View file

@ -106,6 +106,8 @@
<value>log</value>
<value>logarithm</value>
<value>logarithmic</value>
<value>material-union</value>
<value>union</value>
</choice>
</element>
</optional>

8
tests/test_union_energy_grids/geometry.xml Normal file
View file

@ -0,0 +1,8 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" surfaces="-1" />
</geometry>

11
tests/test_union_energy_grids/materials.xml Normal file
View file

@ -0,0 +1,11 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" xs="71c" ao="1.0" />
<nuclide name="H-1" xs="71c" ao="0.5" />
<nuclide name="C-Nat" xs="71c" ao="0.5" />
</material>
</materials>

25
tests/test_union_energy_grids/results.py Normal file
View file

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

2
tests/test_union_energy_grids/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
3.215828E-01 2.966835E-03

18
tests/test_union_energy_grids/settings.xml Normal file
View file

@ -0,0 +1,18 @@
<?xml version="1.0"?>
<settings>
<energy_grid>union</energy_grid>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

59
tests/test_union_energy_grids/test_union_energy_grids.py Normal file
View file

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