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

Fix distribcells and refactor initialization quite a bit

Browse source
This commit is contained in:
Paul Romano 2017-08-17 13:19:01 -05:00
parent 96d35a1acb
commit 9cc06eafae
8 changed files with 494 additions and 524 deletions
Download patch file Download diff file Expand all files Collapse all files

10
src/cmfd_input.F90
View file

@ -441,6 +441,16 @@ contains
call filter_dict % add_key(filt % id, i_filt)
end select
! Initialize filters
do i = i_filt_start, i_filt_end
select type (filt => filters(i) % obj)
type is (SurfaceFilter)
! Don't do anything
class default
call filt % initialize()
end select
end do
! Allocate tallies
err = openmc_extend_tallies(3, i_start, i_end)
cmfd_tallies => tallies(i_start:i_end)

5
src/geometry_header.F90
View file

@ -347,8 +347,7 @@ contains
! temperatures to read (which may be different if interpolation is used)
!===============================================================================
subroutine get_temperatures(cells, nuc_temps, sab_temps)
type(Cell), allocatable, intent(in) :: cells(:)
subroutine get_temperatures(nuc_temps, sab_temps)
type(VectorReal), allocatable, intent(out) :: nuc_temps(:)
type(VectorReal), optional, allocatable, intent(out) :: sab_temps(:)
@ -374,7 +373,7 @@ contains
temperature = cells(i) % sqrtkT(1)**2 / K_BOLTZMANN
end if
i_material = material_dict % get_key(cells(i) % material(j))
i_material = cells(i) % material(j)
associate (mat => materials(i_material))
NUC_NAMES_LOOP: do k = 1, size(mat % names)
! Get index in nuc_temps array

197
src/initialize.F90
View file

@ -12,7 +12,7 @@ module initialize
use dict_header, only: DictIntInt, ElemKeyValueII
use set_header, only: SetInt
use error, only: fatal_error, warning
use geometry, only: neighbor_lists, count_instance, calc_offsets, &
use geometry, only: calc_offsets, &
maximum_levels
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
root_universe
@ -26,7 +26,6 @@ module initialize
use output, only: print_version, write_message, print_usage, &
print_plot
use random_lcg, only: initialize_prng
use state_point, only: load_state_point
use string, only: to_str, starts_with, ends_with, str_to_int
use summary, only: write_summary
use tally_header, only: TallyObject
@ -96,16 +95,9 @@ contains
! XML files because we need the PRNG to be initialized first
if (run_mode == MODE_PLOTTING) call read_plots_xml()
! Use dictionaries to redefine index pointers
call adjust_indices()
! Initialize distribcell_filters
call prepare_distribcell()
! After reading input and basic geometry setup is complete, build lists of
! neighboring cells for efficient tracking
call neighbor_lists()
! Check to make sure there are not too many nested coordinate levels in the
! geometry since the coordinate list is statically allocated for performance
! reasons
@ -126,9 +118,6 @@ contains
! fission bank
call allocate_banks()
! If this is a restart run, load the state point data and binary source
! file
if (restart_run) call load_state_point()
end if
if (master) then
@ -420,176 +409,6 @@ contains
end subroutine read_command_line
!===============================================================================
! ADJUST_INDICES changes the values for 'surfaces' for each cell and the
! material index assigned to each to the indices in the surfaces and material
! array rather than the unique IDs assigned to each surface and material. Also
! assigns boundary conditions to surfaces based on those read into the bc_dict
! dictionary
!===============================================================================
subroutine adjust_indices()
integer :: i ! index for various purposes
integer :: j ! index for various purposes
integer :: k ! loop index for lattices
integer :: m ! loop index for lattices
integer :: lid ! lattice IDs
integer :: i_array ! index in surfaces/materials array
integer :: id ! user-specified id
type(Cell), pointer :: c => null()
class(Lattice), pointer :: lat => null()
do i = 1, n_cells
! =======================================================================
! ADJUST REGION SPECIFICATION FOR EACH CELL
c => cells(i)
do j = 1, size(c%region)
id = c%region(j)
! Make sure that only regions are checked. Since OP_UNION is the
! operator with the lowest integer value, anything below it must denote
! a half-space
if (id < OP_UNION) then
if (surface_dict%has_key(abs(id))) then
i_array = surface_dict%get_key(abs(id))
c%region(j) = sign(i_array, id)
else
call fatal_error("Could not find surface " // trim(to_str(abs(id)))&
&// " specified on cell " // trim(to_str(c%id)))
end if
end if
end do
! Also adjust the indices in the reverse Polish notation
do j = 1, size(c%rpn)
id = c%rpn(j)
! Again, make sure that only regions are checked
if (id < OP_UNION) then
i_array = surface_dict%get_key(abs(id))
c%rpn(j) = sign(i_array, id)
end if
end do
! =======================================================================
! ADJUST UNIVERSE INDEX FOR EACH CELL
id = c%universe
if (universe_dict%has_key(id)) then
c%universe = universe_dict%get_key(id)
else
call fatal_error("Could not find universe " // trim(to_str(id)) &
&// " specified on cell " // trim(to_str(c%id)))
end if
! =======================================================================
! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL
if (c % material(1) == NONE) then
id = c % fill
if (universe_dict % has_key(id)) then
c % type = FILL_UNIVERSE
c % fill = universe_dict % get_key(id)
elseif (lattice_dict % has_key(id)) then
lid = lattice_dict % get_key(id)
c % type = FILL_LATTICE
c % fill = lid
else
call fatal_error("Specified fill " // trim(to_str(id)) // " on cell "&
// trim(to_str(c % id)) // " is neither a universe nor a &
&lattice.")
end if
else
do j = 1, size(c % material)
id = c % material(j)
if (id == MATERIAL_VOID) then
c % type = FILL_MATERIAL
else if (material_dict % has_key(id)) then
c % type = FILL_MATERIAL
c % material(j) = material_dict % get_key(id)
else
call fatal_error("Could not find material " // trim(to_str(id)) &
// " specified on cell " // trim(to_str(c % id)))
end if
end do
end if
end do
! ==========================================================================
! ADJUST UNIVERSE INDICES FOR EACH LATTICE
do i = 1, n_lattices
lat => lattices(i)%obj
select type (lat)
type is (RectLattice)
do m = 1, lat%n_cells(3)
do k = 1, lat%n_cells(2)
do j = 1, lat%n_cells(1)
id = lat%universes(j,k,m)
if (universe_dict%has_key(id)) then
lat%universes(j,k,m) = universe_dict%get_key(id)
else
call fatal_error("Invalid universe number " &
&// trim(to_str(id)) // " specified on lattice " &
&// trim(to_str(lat%id)))
end if
end do
end do
end do
type is (HexLattice)
do m = 1, lat%n_axial
do k = 1, 2*lat%n_rings - 1
do j = 1, 2*lat%n_rings - 1
if (j + k < lat%n_rings + 1) then
cycle
else if (j + k > 3*lat%n_rings - 1) then
cycle
end if
id = lat%universes(j, k, m)
if (universe_dict%has_key(id)) then
lat%universes(j, k, m) = universe_dict%get_key(id)
else
call fatal_error("Invalid universe number " &
&// trim(to_str(id)) // " specified on lattice " &
&// trim(to_str(lat%id)))
end if
end do
end do
end do
end select
if (lat%outer /= NO_OUTER_UNIVERSE) then
if (universe_dict%has_key(lat%outer)) then
lat%outer = universe_dict%get_key(lat%outer)
else
call fatal_error("Invalid universe number " &
&// trim(to_str(lat%outer)) &
&// " specified on lattice " // trim(to_str(lat%id)))
end if
end if
end do
! =======================================================================
! ADJUST INDICES FOR EACH TALLY FILTER
FILTER_LOOP: do i = 1, n_filters
select type(filt => filters(i) % obj)
type is (SurfaceFilter)
! Check if this is a surface filter only for surface currents
if (.not. filt % current) call filt % initialize()
class default
call filt % initialize()
end select
end do FILTER_LOOP
end subroutine adjust_indices
!===============================================================================
! CALCULATE_WORK determines how many particles each processor should simulate
!===============================================================================
@ -717,20 +536,6 @@ contains
! If distribcell isn't used in this simulation then no more work left to do.
if (.not. distribcell_active) return
! Count the number of instances of each cell.
call count_instance(universes(root_universe))
! Set the number of bins in all distribcell filters.
do i = 1, n_tallies
do j = 1, size(tallies(i) % obj % filter)
select type(filt => filters(tallies(i) % obj % filter(j)) % obj)
type is (DistribcellFilter)
! Set the number of bins to the number of instances of the cell.
filt % n_bins = cells(filt % cell) % instances
end select
end do
end do
! Make sure the number of materials and temperatures matches the number of
! cell instances.
do i = 1, n_cells

789
src/input_xml.F90
View file

@ -11,6 +11,7 @@ module input_xml
use distribution_univariate
use endf, only: reaction_name
use error, only: fatal_error, warning
use geometry, only: count_instance, neighbor_lists
use geometry_header, only: Cell, Lattice, RectLattice, HexLattice, &
get_temperatures, root_universe
use global
@ -47,9 +48,25 @@ contains
subroutine read_input_xml()
type(VectorReal), allocatable :: nuc_temps(:) ! List of T to read for each nuclide
type(VectorReal), allocatable :: sab_temps(:) ! List of T to read for each S(a,b)
real(8), allocatable :: material_temps(:)
call read_settings_xml()
call read_cross_sections_xml()
call read_materials_xml(material_temps)
call read_geometry_xml()
call read_materials()
! Set up neighbor lists, convert user IDs -> indices, assign temperatures
call finalize_geometry(material_temps, nuc_temps, sab_temps)
! Read continuous-energy cross sections
if (run_CE .and. run_mode /= MODE_PLOTTING) then
call time_read_xs % start()
call read_ce_cross_sections(nuc_temps, sab_temps)
call time_read_xs % stop()
end if
call read_tallies_xml()
if (cmfd_run) call configure_cmfd()
@ -66,6 +83,27 @@ contains
end subroutine read_input_xml
subroutine finalize_geometry(material_temps, nuc_temps, sab_temps)
real(8), intent(in) :: material_temps(:)
type(VectorReal), allocatable, intent(out) :: nuc_temps(:)
type(VectorReal), optional, allocatable, intent(out) :: sab_temps(:)
! Perform some final operations to set up the geometry
call adjust_indices()
call count_instance(universes(root_universe))
! After reading input and basic geometry setup is complete, build lists of
! neighboring cells for efficient tracking
call neighbor_lists()
! Assign temperatures to cells that don't have temperatures already assigned
call assign_temperatures(material_temps)
! Determine desired txemperatures for each nuclide and S(a,b) table
call get_temperatures(nuc_temps, sab_temps)
end subroutine finalize_geometry
!===============================================================================
! READ_SETTINGS_XML reads data from a settings.xml file and parses it, checking
! for errors and placing properly-formatted data in the right data structures
@ -1078,6 +1116,7 @@ contains
integer :: i, j, k, m, i_x, i_a, input_index
integer :: n, n_mats, n_x, n_y, n_z, n_rings, n_rlats, n_hlats
integer :: id
integer :: univ_id
integer :: n_cells_in_univ
integer :: coeffs_reqd
@ -1111,13 +1150,9 @@ contains
type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
! universe contains
! Display output message
call write_message("Reading geometry XML file...", 5)
! ==========================================================================
! READ CELLS FROM GEOMETRY.XML
! Check if geometry.xml exists
filename = trim(path_input) // "geometry.xml"
inquire(FILE=filename, EXIST=file_exists)
@ -1130,265 +1165,6 @@ contains
call doc % load_file(filename)
root = doc % document_element()
! Get pointer to list of XML <cell>
call get_node_list(root, "cell", node_cell_list)
! Get number of <cell> tags
n_cells = size(node_cell_list)
! Check for no cells
if (n_cells == 0) then
call fatal_error("No cells found in geometry.xml!")
end if
! Allocate cells array
allocate(cells(n_cells))
if (check_overlaps) then
allocate(overlap_check_cnt(n_cells))
overlap_check_cnt = 0
end if
n_universes = 0
do i = 1, n_cells
c => cells(i)
! Initialize distribcell instances and distribcell index
c % instances = 0
c % distribcell_index = NONE
! Get pointer to i-th cell node
node_cell = node_cell_list(i)
! Copy data into cells
if (check_for_node(node_cell, "id")) then
call get_node_value(node_cell, "id", c % id)
else
call fatal_error("Must specify id of cell in geometry XML file.")
end if
! Copy cell name
if (check_for_node(node_cell, "name")) then
call get_node_value(node_cell, "name", c % name)
end if
if (check_for_node(node_cell, "universe")) then
call get_node_value(node_cell, "universe", c % universe)
else
c % universe = 0
end if
if (check_for_node(node_cell, "fill")) then
call get_node_value(node_cell, "fill", c % fill)
if (find(fill_univ_ids, c % fill) == -1) &
call fill_univ_ids % push_back(c % fill)
else
c % fill = NONE
end if
! Check to make sure 'id' hasn't been used
if (cell_dict % has_key(c % id)) then
call fatal_error("Two or more cells use the same unique ID: " &
// to_str(c % id))
end if
! Read material
if (check_for_node(node_cell, "material")) then
n_mats = node_word_count(node_cell, "material")
if (n_mats > 0) then
allocate(sarray(n_mats))
call get_node_array(node_cell, "material", sarray)
allocate(c % material(n_mats))
do j = 1, n_mats
select case(trim(to_lower(sarray(j))))
case ('void')
c % material(j) = MATERIAL_VOID
case default
c % material(j) = int(str_to_int(sarray(j)), 4)
! Check for error
if (c % material(j) == ERROR_INT) then
call fatal_error("Invalid material specified on cell " &
// to_str(c % id))
end if
end select
end do
deallocate(sarray)
else
allocate(c % material(1))
c % material(1) = NONE
end if
else
allocate(c % material(1))
c % material(1) = NONE
end if
! Check to make sure that either material or fill was specified
if (c % material(1) == NONE .and. c % fill == NONE) then
call fatal_error("Neither material nor fill was specified for cell " &
// trim(to_str(c % id)))
end if
! Check to make sure that both material and fill haven't been
! specified simultaneously
if (c % material(1) /= NONE .and. c % fill /= NONE) then
call fatal_error("Cannot specify material and fill simultaneously")
end if
! Check for region specification (also under deprecated name surfaces)
if (check_for_node(node_cell, "surfaces")) then
call warning("The use of 'surfaces' is deprecated and will be &
&disallowed in a future release. Use 'region' instead. The &
&openmc-update-inputs utility can be used to automatically &
&update geometry.xml files.")
region_spec = node_value_string(node_cell, "surfaces")
call get_node_value(node_cell, "surfaces", region_spec)
elseif (check_for_node(node_cell, "region")) then
region_spec = node_value_string(node_cell, "region")
else
region_spec = ''
end if
if (len_trim(region_spec) > 0) then
! Create surfaces array from string
call tokenize(region_spec, tokens)
! Use shunting-yard algorithm to determine RPN for surface algorithm
call generate_rpn(c%id, tokens, rpn)
! Copy region spec and RPN form to cell arrays
allocate(c % region(tokens%size()))
allocate(c % rpn(rpn%size()))
c % region(:) = tokens%data(1:tokens%size())
c % rpn(:) = rpn%data(1:rpn%size())
call tokens%clear()
call rpn%clear()
end if
if (.not. allocated(c%region)) allocate(c%region(0))
if (.not. allocated(c%rpn)) allocate(c%rpn(0))
! Check if this is a simple cell
if (any(c%rpn == OP_COMPLEMENT) .or. any(c%rpn == OP_UNION)) then
c%simple = .false.
else
c%simple = .true.
end if
! Rotation matrix
if (check_for_node(node_cell, "rotation")) then
! Rotations can only be applied to cells that are being filled with
! another universe
if (c % fill == NONE) then
call fatal_error("Cannot apply a rotation to cell " // trim(to_str(&
&c % id)) // " because it is not filled with another universe")
end if
! Read number of rotation parameters
n = node_word_count(node_cell, "rotation")
if (n /= 3) then
call fatal_error("Incorrect number of rotation parameters on cell " &
// to_str(c % id))
end if
! Copy rotation angles in x,y,z directions
allocate(c % rotation(3))
call get_node_array(node_cell, "rotation", c % rotation)
phi = -c % rotation(1) * PI/180.0_8
theta = -c % rotation(2) * PI/180.0_8
psi = -c % rotation(3) * PI/180.0_8
! Calculate rotation matrix based on angles given
allocate(c % rotation_matrix(3,3))
c % rotation_matrix = reshape((/ &
cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), &
-cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), &
cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), &
sin(phi)*cos(theta), &
sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), &
-sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), &
cos(phi)*cos(theta) /), (/ 3,3 /))
end if
! Translation vector
if (check_for_node(node_cell, "translation")) then
! Translations can only be applied to cells that are being filled with
! another universe
if (c % fill == NONE) then
call fatal_error("Cannot apply a translation to cell " &
// trim(to_str(c % id)) // " because it is not filled with &
&another universe")
end if
! Read number of translation parameters
n = node_word_count(node_cell, "translation")
if (n /= 3) then
call fatal_error("Incorrect number of translation parameters on &
&cell " // to_str(c % id))
end if
! Copy translation vector
allocate(c % translation(3))
call get_node_array(node_cell, "translation", c % translation)
end if
! Read cell temperatures. If the temperature is not specified, set it to
! ERROR_REAL for now. During initialization we'll replace ERROR_REAL with
! the temperature from the material data.
if (check_for_node(node_cell, "temperature")) then
n = node_word_count(node_cell, "temperature")
if (n > 0) then
! Make sure this is a "normal" cell.
if (c % material(1) == NONE) call fatal_error("Cell " &
// trim(to_str(c % id)) // " was specified with a temperature &
&but no material. Temperature specification is only valid for &
&cells filled with a material.")
! Copy in temperatures
allocate(c % sqrtkT(n))
call get_node_array(node_cell, "temperature", c % sqrtkT)
! Make sure all temperatues are positive
do j = 1, size(c % sqrtkT)
if (c % sqrtkT(j) < ZERO) call fatal_error("Cell " &
// trim(to_str(c % id)) // " was specified with a negative &
&temperature. All cell temperatures must be non-negative.")
end do
! Convert to sqrt(kT)
c % sqrtkT(:) = sqrt(K_BOLTZMANN * c % sqrtkT(:))
else
allocate(c % sqrtkT(1))
c % sqrtkT(1) = ERROR_REAL
end if
else
allocate(c % sqrtkT(1))
c % sqrtkT = ERROR_REAL
end if
! Add cell to dictionary
call cell_dict % add_key(c % id, i)
! For cells, we also need to check if there's a new universe --
! also for every cell add 1 to the count of cells for the
! specified universe
univ_id = c % universe
if (.not. cells_in_univ_dict % has_key(univ_id)) then
n_universes = n_universes + 1
n_cells_in_univ = 1
call universe_dict % add_key(univ_id, n_universes)
call univ_ids % push_back(univ_id)
else
n_cells_in_univ = 1 + cells_in_univ_dict % get_key(univ_id)
end if
call cells_in_univ_dict % add_key(univ_id, n_cells_in_univ)
end do
! ==========================================================================
! READ SURFACES FROM GEOMETRY.XML
@ -1688,6 +1464,279 @@ contains
end if
end do
! ==========================================================================
! READ CELLS FROM GEOMETRY.XML
! Get pointer to list of XML <cell>
call get_node_list(root, "cell", node_cell_list)
! Get number of <cell> tags
n_cells = size(node_cell_list)
! Check for no cells
if (n_cells == 0) then
call fatal_error("No cells found in geometry.xml!")
end if
! Allocate cells array
allocate(cells(n_cells))
if (check_overlaps) then
allocate(overlap_check_cnt(n_cells))
overlap_check_cnt = 0
end if
n_universes = 0
do i = 1, n_cells
c => cells(i)
! Initialize distribcell instances and distribcell index
c % instances = 0
c % distribcell_index = NONE
! Get pointer to i-th cell node
node_cell = node_cell_list(i)
! Copy data into cells
if (check_for_node(node_cell, "id")) then
call get_node_value(node_cell, "id", c % id)
else
call fatal_error("Must specify id of cell in geometry XML file.")
end if
! Copy cell name
if (check_for_node(node_cell, "name")) then
call get_node_value(node_cell, "name", c % name)
end if
if (check_for_node(node_cell, "universe")) then
call get_node_value(node_cell, "universe", c % universe)
else
c % universe = 0
end if
if (check_for_node(node_cell, "fill")) then
call get_node_value(node_cell, "fill", c % fill)
if (find(fill_univ_ids, c % fill) == -1) &
call fill_univ_ids % push_back(c % fill)
else
c % fill = NONE
end if
! Check to make sure 'id' hasn't been used
if (cell_dict % has_key(c % id)) then
call fatal_error("Two or more cells use the same unique ID: " &
// to_str(c % id))
end if
! Read material
if (check_for_node(node_cell, "material")) then
n_mats = node_word_count(node_cell, "material")
if (n_mats > 0) then
allocate(sarray(n_mats))
call get_node_array(node_cell, "material", sarray)
allocate(c % material(n_mats))
do j = 1, n_mats
select case(trim(to_lower(sarray(j))))
case ('void')
c % material(j) = MATERIAL_VOID
case default
c % material(j) = int(str_to_int(sarray(j)), 4)
! Check for error
if (c % material(j) == ERROR_INT) then
call fatal_error("Invalid material specified on cell " &
// to_str(c % id))
end if
end select
end do
deallocate(sarray)
else
allocate(c % material(1))
c % material(1) = NONE
end if
else
allocate(c % material(1))
c % material(1) = NONE
end if
! Check to make sure that either material or fill was specified
if (c % material(1) == NONE .and. c % fill == NONE) then
call fatal_error("Neither material nor fill was specified for cell " &
// trim(to_str(c % id)))
end if
! Check to make sure that both material and fill haven't been
! specified simultaneously
if (c % material(1) /= NONE .and. c % fill /= NONE) then
call fatal_error("Cannot specify material and fill simultaneously")
end if
! Check for region specification (also under deprecated name surfaces)
if (check_for_node(node_cell, "surfaces")) then
call warning("The use of 'surfaces' is deprecated and will be &
&disallowed in a future release. Use 'region' instead. The &
&openmc-update-inputs utility can be used to automatically &
&update geometry.xml files.")
region_spec = node_value_string(node_cell, "surfaces")
call get_node_value(node_cell, "surfaces", region_spec)
elseif (check_for_node(node_cell, "region")) then
region_spec = node_value_string(node_cell, "region")
else
region_spec = ''
end if
if (len_trim(region_spec) > 0) then
! Create surfaces array from string
call tokenize(region_spec, tokens)
! Convert user IDs to surface indices
do j = 1, tokens % size()
id = tokens % data(j)
if (id < OP_UNION) then
if (surface_dict % has_key(abs(id))) then
k = surface_dict % get_key(abs(id))
tokens % data(j) = sign(k, id)
end if
end if
end do
! Use shunting-yard algorithm to determine RPN for surface algorithm
call generate_rpn(c%id, tokens, rpn)
! Copy region spec and RPN form to cell arrays
allocate(c % region(tokens%size()))
allocate(c % rpn(rpn%size()))
c % region(:) = tokens%data(1:tokens%size())
c % rpn(:) = rpn%data(1:rpn%size())
call tokens%clear()
call rpn%clear()
end if
if (.not. allocated(c%region)) allocate(c%region(0))
if (.not. allocated(c%rpn)) allocate(c%rpn(0))
! Check if this is a simple cell
if (any(c%rpn == OP_COMPLEMENT) .or. any(c%rpn == OP_UNION)) then
c%simple = .false.
else
c%simple = .true.
end if
! Rotation matrix
if (check_for_node(node_cell, "rotation")) then
! Rotations can only be applied to cells that are being filled with
! another universe
if (c % fill == NONE) then
call fatal_error("Cannot apply a rotation to cell " // trim(to_str(&
&c % id)) // " because it is not filled with another universe")
end if
! Read number of rotation parameters
n = node_word_count(node_cell, "rotation")
if (n /= 3) then
call fatal_error("Incorrect number of rotation parameters on cell " &
// to_str(c % id))
end if
! Copy rotation angles in x,y,z directions
allocate(c % rotation(3))
call get_node_array(node_cell, "rotation", c % rotation)
phi = -c % rotation(1) * PI/180.0_8
theta = -c % rotation(2) * PI/180.0_8
psi = -c % rotation(3) * PI/180.0_8
! Calculate rotation matrix based on angles given
allocate(c % rotation_matrix(3,3))
c % rotation_matrix = reshape((/ &
cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), &
-cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), &
cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), &
sin(phi)*cos(theta), &
sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), &
-sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), &
cos(phi)*cos(theta) /), (/ 3,3 /))
end if
! Translation vector
if (check_for_node(node_cell, "translation")) then
! Translations can only be applied to cells that are being filled with
! another universe
if (c % fill == NONE) then
call fatal_error("Cannot apply a translation to cell " &
// trim(to_str(c % id)) // " because it is not filled with &
&another universe")
end if
! Read number of translation parameters
n = node_word_count(node_cell, "translation")
if (n /= 3) then
call fatal_error("Incorrect number of translation parameters on &
&cell " // to_str(c % id))
end if
! Copy translation vector
allocate(c % translation(3))
call get_node_array(node_cell, "translation", c % translation)
end if
! Read cell temperatures. If the temperature is not specified, set it to
! ERROR_REAL for now. During initialization we'll replace ERROR_REAL with
! the temperature from the material data.
if (check_for_node(node_cell, "temperature")) then
n = node_word_count(node_cell, "temperature")
if (n > 0) then
! Make sure this is a "normal" cell.
if (c % material(1) == NONE) call fatal_error("Cell " &
// trim(to_str(c % id)) // " was specified with a temperature &
&but no material. Temperature specification is only valid for &
&cells filled with a material.")
! Copy in temperatures
allocate(c % sqrtkT(n))
call get_node_array(node_cell, "temperature", c % sqrtkT)
! Make sure all temperatues are positive
do j = 1, size(c % sqrtkT)
if (c % sqrtkT(j) < ZERO) call fatal_error("Cell " &
// trim(to_str(c % id)) // " was specified with a negative &
&temperature. All cell temperatures must be non-negative.")
end do
! Convert to sqrt(kT)
c % sqrtkT(:) = sqrt(K_BOLTZMANN * c % sqrtkT(:))
else
allocate(c % sqrtkT(1))
c % sqrtkT(1) = ERROR_REAL
end if
else
allocate(c % sqrtkT(1))
c % sqrtkT = ERROR_REAL
end if
! Add cell to dictionary
call cell_dict % add_key(c % id, i)
! For cells, we also need to check if there's a new universe --
! also for every cell add 1 to the count of cells for the
! specified universe
univ_id = c % universe
if (.not. cells_in_univ_dict % has_key(univ_id)) then
n_universes = n_universes + 1
n_cells_in_univ = 1
call universe_dict % add_key(univ_id, n_universes)
call univ_ids % push_back(univ_id)
else
n_cells_in_univ = 1 + cells_in_univ_dict % get_key(univ_id)
end if
call cells_in_univ_dict % add_key(univ_id, n_cells_in_univ)
end do
! ==========================================================================
! READ LATTICES FROM GEOMETRY.XML
@ -2071,11 +2120,8 @@ contains
! for errors and placing properly-formatted data in the right data structures
!===============================================================================
subroutine read_materials()
subroutine read_cross_sections_xml()
integer :: i, j
type(VectorReal), allocatable :: nuc_temps(:) ! List of T to read for each nuclide
type(VectorReal), allocatable :: sab_temps(:) ! List of T to read for each S(a,b)
real(8), allocatable :: material_temps(:)
logical :: file_exists
character(MAX_FILE_LEN) :: env_variable
character(MAX_LINE_LEN) :: filename
@ -2183,22 +2229,7 @@ contains
end do
end if
! Parse data from materials.xml
call read_materials_xml(material_temps)
! Assign temperatures to cells that don't have temperatures already assigned
call assign_temperatures(material_temps)
! Determine desired temperatures for each nuclide and S(a,b) table
call get_temperatures(cells, nuc_temps, sab_temps)
! Read continuous-energy cross sections
if (run_CE .and. run_mode /= MODE_PLOTTING) then
call time_read_xs % start()
call read_ce_cross_sections(nuc_temps, sab_temps)
call time_read_xs % stop()
end if
end subroutine read_materials
end subroutine read_cross_sections_xml
subroutine read_materials_xml(material_temps)
real(8), allocatable, intent(out) :: material_temps(:)
@ -3315,6 +3346,9 @@ contains
! Set filter id
err = openmc_filter_set_id(i_start + i - 1, filter_id)
! Initialize filter
call f % obj % initialize()
end do READ_FILTERS
! ==========================================================================
@ -3414,6 +3448,7 @@ contains
else
allocate(t % filter(n_filter))
end if
deallocate(temp_filter)
! =======================================================================
! READ DATA FOR NUCLIDES
@ -3822,51 +3857,28 @@ contains
! Get pointer to mesh
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
i_mesh = filt % mesh
m => meshes(filt % mesh)
end select
! Copy filter indices to temporary array
allocate(temp_filter(size(t % filter) + 1))
temp_filter(1:size(t % filter)) = t % filter
! Move allocation back -- temp_filter becomes deallocated during
! this call
call move_alloc(FROM=temp_filter, TO=t % filter)
n_filter = size(t % filter)
! Extend the filters array so we can add a surface filter and
! mesh filter
err = openmc_extend_filters(2, i_filt_start, i_filt_end)
! Get index of the new mesh filter
i_filt = i_filt_start
! Duplicate the mesh filter since other tallies might use this
! filter and we need to change the dimension
allocate(MeshFilter :: filters(i_filt) % obj)
select type(filt => filters(i_filt) % obj)
type is (MeshFilter)
filt % id = i_filt
filt % mesh = i_mesh
filters(i_filt_start) = filters(i_filter_mesh)
! We need to increase the dimension by one since we also need
! currents coming into and out of the boundary mesh cells.
filt % n_bins = product(m % dimension + 1)
! We need to increase the dimension by one since we also need
! currents coming into and out of the boundary mesh cells.
filters(i_filt_start) % obj % n_bins = product(m % dimension + 1)
! Add filter to dictionary
call filter_dict % add_key(filt % id, i_filt)
end select
t % filter(t % find_filter(FILTER_MESH)) = i_filt
! Get index of the new surface filter
i_filt = i_filt_end
! Set ID
err = openmc_filter_set_id(i_filt_start, i_filt_start)
! Add surface filter
allocate(SurfaceFilter :: filters(i_filt) % obj)
select type (filt => filters(i_filt) % obj)
allocate(SurfaceFilter :: filters(i_filt_end) % obj)
select type (filt => filters(i_filt_end) % obj)
type is (SurfaceFilter)
filt % id = i_filt
filt % n_bins = 4 * m % n_dimension
allocate(filt % surfaces(4 * m % n_dimension))
if (m % n_dimension == 1) then
@ -3881,11 +3893,23 @@ contains
end if
filt % current = .true.
! Add filter to dictionary
call filter_dict % add_key(filt % id, i_filt)
! Set ID
err = openmc_filter_set_id(i_filt_end, i_filt_end)
end select
t % find_filter(FILTER_SURFACE) = n_filter
t % filter(n_filter) = i_filt
! Copy filter indices to resized array
n_filter = size(t % filter)
allocate(temp_filter(n_filter + 1))
temp_filter(1:size(t % filter)) = t % filter
n_filter = n_filter + 1
! Set mesh and surface filters
temp_filter(t % find_filter(FILTER_MESH)) = i_filt_start
temp_filter(n_filter) = i_filt_end
! Set filters
err = openmc_tally_set_filters(i_start + i - 1, n_filter, temp_filter)
deallocate(temp_filter)
end if
case ('events')
@ -5289,7 +5313,7 @@ contains
end if
! Use material default or global default temperature
i_material = material_dict % get_key(cells(i) % material(j))
i_material = cells(i) % material(j)
if (material_temps(i_material) /= ERROR_REAL) then
cells(i) % sqrtkT(j) = sqrt(K_BOLTZMANN * &
material_temps(i_material))
@ -5428,4 +5452,129 @@ contains
end if
end subroutine check_data_version
!===============================================================================
! ADJUST_INDICES changes the values for 'surfaces' for each cell and the
! material index assigned to each to the indices in the surfaces and material
! array rather than the unique IDs assigned to each surface and material. Also
! assigns boundary conditions to surfaces based on those read into the bc_dict
! dictionary
!===============================================================================
subroutine adjust_indices()
integer :: i ! index for various purposes
integer :: j ! index for various purposes
integer :: k ! loop index for lattices
integer :: m ! loop index for lattices
integer :: lid ! lattice IDs
integer :: id ! user-specified id
class(Lattice), pointer :: lat => null()
do i = 1, n_cells
! =======================================================================
! ADJUST UNIVERSE INDEX FOR EACH CELL
associate (c => cells(i))
id = c%universe
if (universe_dict%has_key(id)) then
c%universe = universe_dict%get_key(id)
else
call fatal_error("Could not find universe " // trim(to_str(id)) &
&// " specified on cell " // trim(to_str(c%id)))
end if
! =======================================================================
! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL
if (c % material(1) == NONE) then
id = c % fill
if (universe_dict % has_key(id)) then
c % type = FILL_UNIVERSE
c % fill = universe_dict % get_key(id)
elseif (lattice_dict % has_key(id)) then
lid = lattice_dict % get_key(id)
c % type = FILL_LATTICE
c % fill = lid
else
call fatal_error("Specified fill " // trim(to_str(id)) // " on cell "&
// trim(to_str(c % id)) // " is neither a universe nor a &
&lattice.")
end if
else
do j = 1, size(c % material)
id = c % material(j)
if (id == MATERIAL_VOID) then
c % type = FILL_MATERIAL
else if (material_dict % has_key(id)) then
c % type = FILL_MATERIAL
c % material(j) = material_dict % get_key(id)
else
call fatal_error("Could not find material " // trim(to_str(id)) &
// " specified on cell " // trim(to_str(c % id)))
end if
end do
end if
end associate
end do
! ==========================================================================
! ADJUST UNIVERSE INDICES FOR EACH LATTICE
do i = 1, n_lattices
lat => lattices(i)%obj
select type (lat)
type is (RectLattice)
do m = 1, lat%n_cells(3)
do k = 1, lat%n_cells(2)
do j = 1, lat%n_cells(1)
id = lat%universes(j,k,m)
if (universe_dict%has_key(id)) then
lat%universes(j,k,m) = universe_dict%get_key(id)
else
call fatal_error("Invalid universe number " &
&// trim(to_str(id)) // " specified on lattice " &
&// trim(to_str(lat%id)))
end if
end do
end do
end do
type is (HexLattice)
do m = 1, lat%n_axial
do k = 1, 2*lat%n_rings - 1
do j = 1, 2*lat%n_rings - 1
if (j + k < lat%n_rings + 1) then
cycle
else if (j + k > 3*lat%n_rings - 1) then
cycle
end if
id = lat%universes(j, k, m)
if (universe_dict%has_key(id)) then
lat%universes(j, k, m) = universe_dict%get_key(id)
else
call fatal_error("Invalid universe number " &
&// trim(to_str(id)) // " specified on lattice " &
&// trim(to_str(lat%id)))
end if
end do
end do
end do
end select
if (lat%outer /= NO_OUTER_UNIVERSE) then
if (universe_dict%has_key(lat%outer)) then
lat%outer = universe_dict%get_key(lat%outer)
else
call fatal_error("Invalid universe number " &
&// trim(to_str(lat%outer)) &
&// " specified on lattice " // trim(to_str(lat%id)))
end if
end if
end do
end subroutine adjust_indices
end module input_xml

4
src/mgxs_data.F90
View file

@ -51,7 +51,7 @@ contains
call write_message("Loading cross section data...", 5)
! Get temperatures
call get_temperatures(cells, temps)
call get_temperatures(temps)
! Open file for reading
file_id = file_open(path_cross_sections, 'r', parallel=.true.)
@ -233,7 +233,7 @@ contains
kT = cells(i) % sqrtkT(1)**2
end if
i_material = material_dict % get_key(cells(i) % material(j))
i_material = cells(i) % material(j)
! Add temperature if it hasn't already been added
if (find(kTs(i_material), kT) == -1) then

10
src/simulation.F90
View file

@ -18,7 +18,7 @@ module simulation
use particle_header, only: Particle
use random_lcg, only: set_particle_seed
use source, only: initialize_source, sample_external_source
use state_point, only: write_state_point, write_source_point
use state_point, only: write_state_point, write_source_point, load_state_point
use string, only: to_str
use tally, only: accumulate_tallies, setup_active_tallies
use tally_header, only: configure_tallies
@ -383,7 +383,13 @@ contains
allocate(filter_matches(n_filters))
!$omp end parallel
if (.not. restart_run) call initialize_source()
! If this is a restart run, load the state point data and binary source
! file
if (restart_run) then
call load_state_point()
else
call initialize_source()
end if
! Display header
if (master) then

1
src/tallies/tally_filter_distribcell.F90
View file

@ -85,6 +85,7 @@ contains
id = this % cell
if (cell_dict % has_key(id)) then
this % cell = cell_dict % get_key(id)
this % n_bins = cells(this % cell) % instances
else
call fatal_error("Could not find cell " // trim(to_str(id)) &
&// " specified on tally filter.")

2
src/tallies/tally_header.F90
View file

@ -89,7 +89,7 @@ module tally_header
! second dimension of the array is for the combination of filters
! (e.g. specific cell, specific energy group, etc.)
integer :: n_filter_bins
integer :: n_filter_bins = 1
integer :: total_score_bins
real(C_DOUBLE), allocatable :: results(:,:,:)