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

Refactored rest of code to use new HDF5_interface API

Browse source
This commit is contained in:
Paul Romano 2015-09-03 09:20:47 +07:00
parent 205be92897
commit 2f18e4b57e
10 changed files with 1146 additions and 3008 deletions
Download patch file Download diff file Expand all files Collapse all files

11
src/finalize.F90
View file

@ -9,7 +9,8 @@ module finalize
use message_passing
#endif
use hdf5_interface, only: h5tclose_f, h5close_f, hdf5_err
use hdf5_interface, only: hdf5_bank_t, hdf5_tallyresult_t
use hdf5, only: h5tclose_f, h5close_f
implicit none
@ -22,8 +23,10 @@ contains
subroutine finalize_run()
integer :: hdf5_err
! Start finalization timer
call time_finalize % start()
call time_finalize%start()
if (run_mode /= MODE_PLOTTING .and. run_mode /= MODE_PARTICLE) then
! Calculate statistics for tallies and write to tallies.out
@ -37,8 +40,8 @@ contains
end if
! Stop timers and show timing statistics
call time_finalize % stop()
call time_total % stop()
call time_finalize%stop()
call time_total%stop()
if (master .and. (run_mode /= MODE_PLOTTING .and. &
run_mode /= MODE_PARTICLE)) then
call print_runtime()

222
src/hdf5_interface.F90
View file

@ -64,10 +64,9 @@ module hdf5_interface
public :: file_create
public :: file_open
public :: file_close
public :: create_group
public :: open_group
public :: close_group
public :: write_source_bank
public :: read_source_bank
public :: write_attribute_string
contains
@ -243,7 +242,9 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
@ -287,9 +288,10 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -344,7 +346,9 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
@ -404,9 +408,10 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -461,7 +466,9 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
@ -521,9 +528,10 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -578,7 +586,9 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
@ -638,9 +648,10 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -695,7 +706,9 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
@ -755,9 +768,10 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -795,7 +809,9 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
@ -839,9 +855,10 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -896,7 +913,9 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
@ -956,9 +975,10 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -1013,7 +1033,9 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
@ -1073,9 +1095,10 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -1130,7 +1153,9 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
@ -1190,9 +1215,10 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -1247,7 +1273,9 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
@ -1307,9 +1335,10 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -1347,7 +1376,9 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
@ -1391,9 +1422,10 @@ contains
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -1432,12 +1464,13 @@ contains
integer :: n
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
integer(HSIZE_T) :: dims1(1)
integer(HSIZE_T) :: dims2(2)
type(c_ptr) :: f_ptr
character(len=len_trim(buffer)), dimension(1) :: str_tmp
! Set up collective vs. independent I/O
@ -1492,12 +1525,13 @@ contains
integer :: n
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
integer(HSIZE_T) :: dims1(1)
integer(HSIZE_T) :: dims2(2)
type(c_ptr) :: f_ptr
character(len=len_trim(buffer)), dimension(1) :: str_tmp
! Set up collective vs. independent I/O
@ -1673,164 +1707,6 @@ contains
call h5dclose_f(dset, hdf5_err)
end subroutine read_tally_result_2D_explicit
!===============================================================================
! WRITE_SOURCE_BANK writes OpenMC source_bank data
!===============================================================================
subroutine write_source_bank(group_id)
use bank_header, only: Bank
use global, only: n_particles, work, source_bank
integer(HID_T), intent(in) :: group_id
integer :: hdf5_err
integer :: data_xfer_mode
integer(HID_T) :: plist ! property list
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
integer(HID_T) :: memspace ! memory space handle
integer(HSIZE_T) :: dims(1)
type(c_ptr) :: f_ptr
#ifdef PHDF5
integer(HSIZE_T) :: offset(1) ! source data offset
#endif
#ifdef PHDF5
! Set size of total dataspace for all procs and rank
dims(1) = n_particles
call h5screate_simple_f(1, dims, dspace, hdf5_err)
call h5dcreate_f(group_id, "source_bank", hdf5_bank_t, dspace, dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
! Create another data space but for each proc individually
dims(1) = work
call h5screate_simple_f(rank, dims, memspace, hdf5_err)
! Get the individual local proc dataspace
call h5dget_space_f(dset, dspace, hdf5_err)
! Select hyperslab for this dataspace
offset(1) = work_index(rank)
call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims, hdf5_err)
! Set up the property list for parallel writing
call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err)
call h5pset_dxpl_mpio_f(plist, H5FD_MPIO_COLLECTIVE_F, hdf5_err)
! Set up pointer to data
f_ptr = c_loc(source_bank)
! Write data to file in parallel
call h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err, &
file_space_id=dspace, mem_space_id=memspace, &
xfer_prp=plist)
! Close all ids
call h5sclose_f(dspace, hdf5_err)
call h5sclose_f(memspace, hdf5_err)
call h5dclose_f(dset, hdf5_err)
call h5pclose_f(plist, hdf5_err)
#else
! Set size
dims(1) = work
! Create dataspace
call h5screate_simple_f(1, dims, dspace, hdf5_err)
! Create dataset
call h5dcreate_f(group_id, "source_bank", hdf5_bank_t, &
dspace, dset, hdf5_err)
! Set up pointer to data
f_ptr = c_loc(source_bank)
! Write dataset to file
call h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err)
! Close all ids
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
#endif
end subroutine write_source_bank
!===============================================================================
! READ_SOURCE_BANK reads OpenMC source_bank data
!===============================================================================
subroutine read_source_bank(group_id)
use bank_header, only: Bank
use global, only: work, source_bank
integer(HID_T), intent(in) :: group_id
integer :: hdf5_err
integer :: data_xfer_mode
integer(HID_T) :: plist ! property list
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data space handle
integer(HID_T) :: memspace ! memory space handle
integer(HSIZE_T) :: dims(1)
type(c_ptr) :: f_ptr
#ifdef PHDF5
integer(HSIZE_T) :: offset(1) ! offset of data
#endif
#ifdef PHDF5
! Open the dataset
call h5dopen_f(group_id, "source_bank", dset, hdf5_err)
! Create another data space but for each proc individually
dims(1) = work
call h5screate_simple_f(1, dims, memspace, hdf5_err)
! Get the individual local proc dataspace
call h5dget_space_f(dset, dspace, hdf5_err)
! Select hyperslab for this dataspace
offset(1) = work_index(rank)
call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims, hdf5_err)
! Set up the property list for parallel writing
call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err)
call h5pset_dxpl_mpio_f(plist, H5FD_MPIO_COLLECTIVE_F, hdf5_err)
! Set up pointer to data
f_ptr = c_loc(source_bank)
! Read data from file in parallel
call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err, &
file_space_id=dspace, mem_space_id=memspace, &
xfer_prp=plist)
! Close all ids
call h5sclose_f(dspace, hdf5_err)
call h5sclose_f(memspace, hdf5_err)
call h5dclose_f(dset, hdf5_err)
call h5pclose_f(plist, hdf5_err)
#else
! Open dataset
call h5dopen_f(group_id, "source_bank", dset, hdf5_err)
! Set up pointer to data
f_ptr = c_loc(source_bank)
! Read dataset from file
call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err)
! Close all ids
call h5dclose_f(dset, hdf5_err)
#endif
end subroutine read_source_bank
function using_mpio_device(obj_id) result(mpio)
integer(HID_T), intent(in) :: obj_id
logical :: mpio

814
src/hdf5_summary.F90
View file

File diff suppressed because it is too large Load diff

313
src/initialize.F90
View file

@ -12,12 +12,14 @@ module initialize
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
&BASE_UNIVERSE
use global
use hdf5_interface, only: file_open, read_dataset, file_close, hdf5_bank_t,&
hdf5_tallyresult_t, hdf5_integer8_t
use hdf5_summary, only: hdf5_write_summary
use input_xml, only: read_input_xml, read_cross_sections_xml, &
cells_in_univ_dict, read_plots_xml
use material_header, only: Material
use output, only: title, header, print_version, write_message, &
print_usage, write_xs_summary, print_plot
use output_interface
use random_lcg, only: initialize_prng
use state_point, only: load_state_point
use string, only: to_str, str_to_int, starts_with, ends_with
@ -32,8 +34,7 @@ module initialize
use omp_lib
#endif
use hdf5_interface
use hdf5_summary, only: hdf5_write_summary
use hdf5
implicit none
@ -49,8 +50,8 @@ contains
subroutine initialize_run()
! Start total and initialization timer
call time_total % start()
call time_initialize % start()
call time_total%start()
call time_initialize%start()
#ifdef MPI
! Setup MPI
@ -110,9 +111,9 @@ contains
call normalize_ao()
! Read ACE-format cross sections
call time_read_xs % start()
call time_read_xs%start()
call read_xs()
call time_read_xs % stop()
call time_read_xs%stop()
! Create linked lists for multiple instances of the same nuclide
call same_nuclide_list()
@ -122,9 +123,9 @@ contains
case (GRID_NUCLIDE)
continue
case (GRID_MAT_UNION)
call time_unionize % start()
call time_unionize%start()
call unionized_grid()
call time_unionize % stop()
call time_unionize%stop()
case (GRID_LOGARITHM)
call logarithmic_grid()
end select
@ -167,7 +168,7 @@ contains
end if
! Stop initialization timer
call time_initialize % stop()
call time_initialize%stop()
end subroutine initialize_run
@ -220,10 +221,10 @@ contains
! CREATE MPI_BANK TYPE
! Determine displacements for MPI_BANK type
call MPI_GET_ADDRESS(b % wgt, bank_disp(1), mpi_err)
call MPI_GET_ADDRESS(b % xyz, bank_disp(2), mpi_err)
call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err)
call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err)
call MPI_GET_ADDRESS(b%wgt, bank_disp(1), mpi_err)
call MPI_GET_ADDRESS(b%xyz, bank_disp(2), mpi_err)
call MPI_GET_ADDRESS(b%uvw, bank_disp(3), mpi_err)
call MPI_GET_ADDRESS(b%E, bank_disp(4), mpi_err)
! Adjust displacements
bank_disp = bank_disp - bank_disp(1)
@ -239,8 +240,8 @@ contains
! CREATE MPI_TALLYRESULT TYPE
! Determine displacements for MPI_BANK type
call MPI_GET_ADDRESS(tr % value, result_base_disp, mpi_err)
call MPI_GET_ADDRESS(tr % sum, result_disp(1), mpi_err)
call MPI_GET_ADDRESS(tr%value, result_base_disp, mpi_err)
call MPI_GET_ADDRESS(tr%sum, result_disp(1), mpi_err)
! Adjust displacements
result_disp = result_disp - result_base_disp
@ -274,6 +275,7 @@ contains
type(TallyResult), target :: tmp(2) ! temporary TallyResult
type(Bank), target :: tmpb(2) ! temporary Bank
integer :: hdf5_err
integer(HID_T) :: coordinates_t ! HDF5 type for 3 reals
integer(HSIZE_T) :: dims(1) = (/3/) ! size of coordinates
@ -318,8 +320,8 @@ contains
integer :: argc ! number of command line arguments
integer :: last_flag ! index of last flag
integer :: filetype
integer(HID_T) :: file_id
character(MAX_WORD_LEN), allocatable :: argv(:) ! command line arguments
type(BinaryOutput) :: sp
! Check number of command line arguments and allocate argv
argc = COMMAND_ARGUMENT_COUNT()
@ -356,9 +358,9 @@ contains
i = i + 1
! Check what type of file this is
call sp % file_open(argv(i), 'r', serial = .false.)
call sp % read_data(filetype, 'filetype')
call sp % file_close()
file_id = file_open(argv(i), 'r', parallel=.true.)
call read_dataset(file_id, 'filetype', filetype)
call file_close(file_id)
! Set path and flag for type of run
select case (filetype)
@ -379,13 +381,12 @@ contains
i = i + 1
! Check if it has extension we can read
if ((ends_with(argv(i), '.binary') .or. &
ends_with(argv(i), '.h5'))) then
if (ends_with(argv(i), '.h5')) then
! Check file type is a source file
call sp % file_open(argv(i), 'r', serial = .false.)
call sp % read_data(filetype, 'filetype')
call sp % file_close()
file_id = file_open(argv(i), 'r', parallel=.true.)
call read_dataset(file_id, 'filetype', filetype)
call file_close(file_id)
if (filetype /= FILETYPE_SOURCE) then
call fatal_error("Second file after restart flag must be a &
&source file")
@ -494,26 +495,26 @@ contains
! pairs are the id of the universe and the index in the array. In
! cells_in_univ_dict, it's the id of the universe and the number of cells.
pair_list => universe_dict % keys()
pair_list => universe_dict%keys()
current => pair_list
do while (associated(current))
! Find index of universe in universes array
i_univ = current % value
i_univ = current%value
univ => universes(i_univ)
univ % id = current % key
univ%id = current%key
! Check for lowest level universe
if (univ % id == 0) BASE_UNIVERSE = i_univ
if (univ%id == 0) BASE_UNIVERSE = i_univ
! Find cell count for this universe
n_cells_in_univ = cells_in_univ_dict % get_key(univ % id)
n_cells_in_univ = cells_in_univ_dict%get_key(univ%id)
! Allocate cell list for universe
allocate(univ % cells(n_cells_in_univ))
univ % n_cells = n_cells_in_univ
allocate(univ%cells(n_cells_in_univ))
univ%n_cells = n_cells_in_univ
! Move to next universe
next => current % next
next => current%next
deallocate(current)
current => next
end do
@ -528,17 +529,17 @@ contains
c => cells(i)
! Get pointer to corresponding universe
i_univ = universe_dict % get_key(c % universe)
i_univ = universe_dict%get_key(c%universe)
univ => universes(i_univ)
! Increment the index for the cells array within the Universe object and
! then store the index of the Cell object in that array
index_cell_in_univ(i_univ) = index_cell_in_univ(i_univ) + 1
univ % cells(index_cell_in_univ(i_univ)) = i
univ%cells(index_cell_in_univ(i_univ)) = i
end do
! Clear dictionary
call cells_in_univ_dict % clear()
call cells_in_univ_dict%clear()
end subroutine prepare_universes
@ -568,15 +569,15 @@ contains
! ADJUST SURFACE LIST FOR EACH CELL
c => cells(i)
do j = 1, c % n_surfaces
id = c % surfaces(j)
do j = 1, c%n_surfaces
id = c%surfaces(j)
if (id < OP_DIFFERENCE) then
if (surface_dict % has_key(abs(id))) then
i_array = surface_dict % get_key(abs(id))
c % surfaces(j) = sign(i_array, id)
if (surface_dict%has_key(abs(id))) then
i_array = surface_dict%get_key(abs(id))
c%surfaces(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)))
&// " specified on cell " // trim(to_str(c%id)))
end if
end if
end do
@ -584,40 +585,40 @@ contains
! =======================================================================
! ADJUST UNIVERSE INDEX FOR EACH CELL
id = c % universe
if (universe_dict % has_key(id)) then
c % universe = universe_dict % get_key(id)
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)))
&// " specified on cell " // trim(to_str(c%id)))
end if
! =======================================================================
! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL
id = c % material
id = c%material
if (id == MATERIAL_VOID) then
c % type = CELL_NORMAL
c%type = CELL_NORMAL
elseif (id /= 0) then
if (material_dict % has_key(id)) then
c % type = CELL_NORMAL
c % material = material_dict % get_key(id)
if (material_dict%has_key(id)) then
c%type = CELL_NORMAL
c%material = 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)))
&// " specified on cell " // trim(to_str(c%id)))
end if
else
id = c % fill
if (universe_dict % has_key(id)) then
c % type = CELL_FILL
c % fill = universe_dict % get_key(id)
elseif (lattice_dict % has_key(id)) then
lid = lattice_dict % get_key(id)
c % type = CELL_LATTICE
c % fill = lid
id = c%fill
if (universe_dict%has_key(id)) then
c%type = CELL_FILL
c%fill = universe_dict%get_key(id)
elseif (lattice_dict%has_key(id)) then
lid = lattice_dict%get_key(id)
c%type = CELL_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 &
&// trim(to_str(c%id)) // " is neither a universe nor a &
&lattice.")
end if
end if
@ -627,41 +628,41 @@ contains
! ADJUST UNIVERSE INDICES FOR EACH LATTICE
do i = 1, n_lattices
lat => lattices(i) % obj
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)
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)))
&// 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
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
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)
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)))
&// trim(to_str(lat%id)))
end if
end do
end do
@ -669,13 +670,13 @@ contains
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)
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)))
&// trim(to_str(lat%outer)) &
&// " specified on lattice " // trim(to_str(lat%id)))
end if
end if
@ -687,68 +688,68 @@ contains
! =======================================================================
! ADJUST INDICES FOR EACH TALLY FILTER
FILTER_LOOP: do j = 1, t % n_filters
FILTER_LOOP: do j = 1, t%n_filters
select case (t % filters(j) % type)
select case (t%filters(j)%type)
case (FILTER_DISTRIBCELL)
do k = 1, size(t % filters(j) % int_bins)
id = t % filters(j) % int_bins(k)
if (cell_dict % has_key(id)) then
t % filters(j) % int_bins(k) = cell_dict % get_key(id)
do k = 1, size(t%filters(j)%int_bins)
id = t%filters(j)%int_bins(k)
if (cell_dict%has_key(id)) then
t%filters(j)%int_bins(k) = cell_dict%get_key(id)
else
call fatal_error("Could not find cell " // trim(to_str(id)) // &
" specified on tally " // trim(to_str(t % id)))
" specified on tally " // trim(to_str(t%id)))
end if
end do
case (FILTER_CELL, FILTER_CELLBORN)
do k = 1, t % filters(j) % n_bins
id = t % filters(j) % int_bins(k)
if (cell_dict % has_key(id)) then
t % filters(j) % int_bins(k) = cell_dict % get_key(id)
do k = 1, t%filters(j)%n_bins
id = t%filters(j)%int_bins(k)
if (cell_dict%has_key(id)) then
t%filters(j)%int_bins(k) = cell_dict%get_key(id)
else
call fatal_error("Could not find cell " // trim(to_str(id)) &
&// " specified on tally " // trim(to_str(t % id)))
&// " specified on tally " // trim(to_str(t%id)))
end if
end do
case (FILTER_SURFACE)
! Check if this is a surface filter only for surface currents
if (any(t % score_bins == SCORE_CURRENT)) cycle FILTER_LOOP
if (any(t%score_bins == SCORE_CURRENT)) cycle FILTER_LOOP
do k = 1, t % filters(j) % n_bins
id = t % filters(j) % int_bins(k)
if (surface_dict % has_key(id)) then
t % filters(j) % int_bins(k) = surface_dict % get_key(id)
do k = 1, t%filters(j)%n_bins
id = t%filters(j)%int_bins(k)
if (surface_dict%has_key(id)) then
t%filters(j)%int_bins(k) = surface_dict%get_key(id)
else
call fatal_error("Could not find surface " // trim(to_str(id)) &
&// " specified on tally " // trim(to_str(t % id)))
&// " specified on tally " // trim(to_str(t%id)))
end if
end do
case (FILTER_UNIVERSE)
do k = 1, t % filters(j) % n_bins
id = t % filters(j) % int_bins(k)
if (universe_dict % has_key(id)) then
t % filters(j) % int_bins(k) = universe_dict % get_key(id)
do k = 1, t%filters(j)%n_bins
id = t%filters(j)%int_bins(k)
if (universe_dict%has_key(id)) then
t%filters(j)%int_bins(k) = universe_dict%get_key(id)
else
call fatal_error("Could not find universe " // trim(to_str(id)) &
&// " specified on tally " // trim(to_str(t % id)))
&// " specified on tally " // trim(to_str(t%id)))
end if
end do
case (FILTER_MATERIAL)
do k = 1, t % filters(j) % n_bins
id = t % filters(j) % int_bins(k)
if (material_dict % has_key(id)) then
t % filters(j) % int_bins(k) = material_dict % get_key(id)
do k = 1, t%filters(j)%n_bins
id = t%filters(j)%int_bins(k)
if (material_dict%has_key(id)) then
t%filters(j)%int_bins(k) = material_dict%get_key(id)
else
call fatal_error("Could not find material " // trim(to_str(id)) &
&// " specified on tally " // trim(to_str(t % id)))
&// " specified on tally " // trim(to_str(t%id)))
end if
end do
@ -786,46 +787,46 @@ contains
do i = 1, n_materials
mat => materials(i)
percent_in_atom = (mat % atom_density(1) > ZERO)
density_in_atom = (mat % density > ZERO)
percent_in_atom = (mat%atom_density(1) > ZERO)
density_in_atom = (mat%density > ZERO)
sum_percent = ZERO
do j = 1, mat % n_nuclides
do j = 1, mat%n_nuclides
! determine atomic weight ratio
index_list = xs_listing_dict % get_key(mat % names(j))
awr = xs_listings(index_list) % awr
index_list = xs_listing_dict%get_key(mat%names(j))
awr = xs_listings(index_list)%awr
! if given weight percent, convert all values so that they are divided
! by awr. thus, when a sum is done over the values, it's actually
! sum(w/awr)
if (.not. percent_in_atom) then
mat % atom_density(j) = -mat % atom_density(j) / awr
mat%atom_density(j) = -mat%atom_density(j) / awr
end if
end do
! determine normalized atom percents. if given atom percents, this is
! straightforward. if given weight percents, the value is w/awr and is
! divided by sum(w/awr)
sum_percent = sum(mat % atom_density)
mat % atom_density = mat % atom_density / sum_percent
sum_percent = sum(mat%atom_density)
mat%atom_density = mat%atom_density / sum_percent
! Change density in g/cm^3 to atom/b-cm. Since all values are now in atom
! percent, the sum needs to be re-evaluated as 1/sum(x*awr)
if (.not. density_in_atom) then
sum_percent = ZERO
do j = 1, mat % n_nuclides
index_list = xs_listing_dict % get_key(mat % names(j))
awr = xs_listings(index_list) % awr
x = mat % atom_density(j)
do j = 1, mat%n_nuclides
index_list = xs_listing_dict%get_key(mat%names(j))
awr = xs_listings(index_list)%awr
x = mat%atom_density(j)
sum_percent = sum_percent + x*awr
end do
sum_percent = ONE / sum_percent
mat % density = -mat % density * N_AVOGADRO &
mat%density = -mat%density * N_AVOGADRO &
/ MASS_NEUTRON * sum_percent
end if
! Calculate nuclide atom densities
mat % atom_density = mat % density * mat % atom_density
mat%atom_density = mat%density * mat%atom_density
end do
end subroutine normalize_ao
@ -942,16 +943,16 @@ contains
! Get pointer to tally
tally => tallies(i)
n_filt = tally % n_filters
n_filt = tally%n_filters
! Loop over the filters to determine how many additional filters
! need to be added to this tally
do j = 1, tally % n_filters
do j = 1, tally%n_filters
! Determine type of filter
if (tally % filters(j) % type == FILTER_DISTRIBCELL) then
if (tally%filters(j)%type == FILTER_DISTRIBCELL) then
count_all = .true.
if (size(tally % filters(j) % int_bins) > 1) then
if (size(tally%filters(j)%int_bins) > 1) then
call fatal_error("A distribcell filter was specified with &
&multiple bins. This feature is not supported.")
end if
@ -975,12 +976,12 @@ contains
tally => tallies(i)
! Initialize the filters
do j = 1, tally % n_filters
do j = 1, tally%n_filters
! Set the number of bins to the number of instances of the cell
if (tally % filters(j) % type == FILTER_DISTRIBCELL) then
c => cells(tally % filters(j) % int_bins(1))
tally % filters(j) % n_bins = c % instances
if (tally%filters(j)%type == FILTER_DISTRIBCELL) then
c => cells(tally%filters(j)%int_bins(1))
tally%filters(j)%n_bins = c%instances
end if
end do
@ -1031,12 +1032,12 @@ contains
do i = 1, n_tallies
tally => tallies(i)
do j = 1, tally % n_filters
filter => tally % filters(j)
do j = 1, tally%n_filters
filter => tally%filters(j)
if (filter % type == FILTER_DISTRIBCELL) then
if (.not. cell_list % contains(filter % int_bins(1))) then
call cell_list % add(filter % int_bins(1))
if (filter%type == FILTER_DISTRIBCELL) then
if (.not. cell_list%contains(filter%int_bins(1))) then
call cell_list%add(filter%int_bins(1))
end if
end if
@ -1047,8 +1048,8 @@ contains
! to determine the number of offset tables to allocate
do i = 1, n_universes
univ => universes(i)
do j = 1, univ % n_cells
if (cell_list % contains(univ % cells(j))) then
do j = 1, univ%n_cells
if (cell_list%contains(univ%cells(j))) then
n_maps = n_maps + 1
end if
end do
@ -1070,29 +1071,29 @@ contains
do i = 1, n_universes
univ => universes(i)
do j = 1, univ % n_cells
do j = 1, univ%n_cells
if (cell_list % contains(univ % cells(j))) then
if (cell_list%contains(univ%cells(j))) then
! Loop over all tallies
do l = 1, n_tallies
tally => tallies(l)
do m = 1, tally % n_filters
filter => tally % filters(m)
do m = 1, tally%n_filters
filter => tally%filters(m)
! Loop over only distribcell filters
! If filter points to cell we just found, set offset index
if (filter % type == FILTER_DISTRIBCELL) then
if (filter % int_bins(1) == univ % cells(j)) then
filter % offset = k
if (filter%type == FILTER_DISTRIBCELL) then
if (filter%int_bins(1) == univ%cells(j)) then
filter%offset = k
end if
end if
end do
end do
univ_list(k) = univ % id
univ_list(k) = univ%id
k = k + 1
end if
end do
@ -1100,26 +1101,26 @@ contains
! Allocate the offset tables for lattices
do i = 1, n_lattices
lat => lattices(i) % obj
lat => lattices(i)%obj
select type(lat)
type is (RectLattice)
allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), &
lat % n_cells(3)))
allocate(lat%offset(n_maps, lat%n_cells(1), lat%n_cells(2), &
lat%n_cells(3)))
type is (HexLattice)
allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, &
2 * lat % n_rings - 1, lat % n_axial))
allocate(lat%offset(n_maps, 2 * lat%n_rings - 1, &
2 * lat%n_rings - 1, lat%n_axial))
end select
lat % offset(:, :, :, :) = 0
lat%offset(:, :, :, :) = 0
end do
! Allocate offset table for fill cells
do i = 1, n_cells
if (cells(i) % material == NONE) then
allocate(cells(i) % offset(n_maps))
if (cells(i)%material == NONE) then
allocate(cells(i)%offset(n_maps))
end if
end do

1725
src/output_interface.F90
View file

File diff suppressed because it is too large Load diff

53
src/particle_restart.F90
View file

@ -6,18 +6,18 @@ module particle_restart
use constants
use geometry_header, only: BASE_UNIVERSE
use global
use hdf5_interface, only: file_open, file_close, read_dataset
use output, only: write_message, print_particle
use output_interface, only: BinaryOutput
use particle_header, only: Particle
use random_lcg, only: set_particle_seed
use tracking, only: transport
use hdf5, only: HID_T
implicit none
private
public :: run_particle_restart
type(BinaryOutput) :: pr ! Binary file
contains
!===============================================================================
@ -34,7 +34,7 @@ contains
verbosity = 10
! Initialize the particle to be tracked
call p % initialize()
call p%initialize()
! Read in the restart information
call read_particle_restart(p, previous_run_mode)
@ -46,9 +46,9 @@ contains
select case (previous_run_mode)
case (MODE_EIGENVALUE)
particle_seed = ((current_batch - 1)*gen_per_batch + &
current_gen - 1)*n_particles + p % id
current_gen - 1)*n_particles + p%id
case (MODE_FIXEDSOURCE)
particle_seed = p % id
particle_seed = p%id
end select
call set_particle_seed(particle_seed)
@ -66,40 +66,41 @@ contains
!===============================================================================
subroutine read_particle_restart(p, previous_run_mode)
type(Particle), intent(inout) :: p
integer, intent(inout) :: previous_run_mode
integer :: int_scalar
integer, intent(inout) :: previous_run_mode
type(Particle), intent(inout) :: p
integer(HID_T) :: file_id
! Write meessage
call write_message("Loading particle restart file " &
&// trim(path_particle_restart) // "...", 1)
! Open file
call pr % file_open(path_particle_restart, 'r')
file_id = file_open(path_particle_restart, 'r')
! Read data from file
call pr % read_data(int_scalar, 'filetype')
call pr % read_data(int_scalar, 'revision')
call pr % read_data(current_batch, 'current_batch')
call pr % read_data(gen_per_batch, 'gen_per_batch')
call pr % read_data(current_gen, 'current_gen')
call pr % read_data(n_particles, 'n_particles')
call pr % read_data(previous_run_mode, 'run_mode')
call pr % read_data(p % id, 'id')
call pr % read_data(p % wgt, 'weight')
call pr % read_data(p % E, 'energy')
call pr % read_data(p % coord(1) % xyz, 'xyz', length=3)
call pr % read_data(p % coord(1) % uvw, 'uvw', length=3)
call read_dataset(file_id, 'filetype', int_scalar)
call read_dataset(file_id, 'revision', int_scalar)
call read_dataset(file_id, 'current_batch', current_batch)
call read_dataset(file_id, 'gen_per_batch', gen_per_batch)
call read_dataset(file_id, 'current_gen', current_gen)
call read_dataset(file_id, 'n_particles', n_particles)
call read_dataset(file_id, 'run_mode', previous_run_mode)
call read_dataset(file_id, 'id', p%id)
call read_dataset(file_id, 'weight', p%wgt)
call read_dataset(file_id, 'energy', p%E)
call read_dataset(file_id, 'xyz', p%coord(1)%xyz)
call read_dataset(file_id, 'uvw', p%coord(1)%uvw)
! Set particle last attributes
p % last_wgt = p % wgt
p % last_xyz = p % coord(1) % xyz
p % last_uvw = p % coord(1) % uvw
p % last_E = p % E
p%last_wgt = p%wgt
p%last_xyz = p%coord(1)%xyz
p%last_uvw = p%coord(1)%uvw
p%last_E = p%E
! Close hdf5 file
call pr % file_close()
call file_close(file_id)
end subroutine read_particle_restart

42
src/particle_restart_write.F90
View file

@ -2,17 +2,16 @@ module particle_restart_write
use bank_header, only: Bank
use global
use output_interface, only: BinaryOutput
use hdf5_interface
use particle_header, only: Particle
use string, only: to_str
use hdf5
implicit none
private
public :: write_particle_restart
! Binary output file
type(BinaryOutput) :: pr
contains
!===============================================================================
@ -20,43 +19,42 @@ contains
!===============================================================================
subroutine write_particle_restart(p)
type(Particle), intent(in) :: p
integer(HID_T) :: file_id
character(MAX_FILE_LEN) :: filename
type(Bank), pointer :: src => null()
type(Bank), pointer :: src
! Dont write another restart file if in particle restart mode
if (run_mode == MODE_PARTICLE) return
! Set up file name
filename = trim(path_output) // 'particle_' // trim(to_str(current_batch)) &
// '_' // trim(to_str(p % id))
filename = trim(filename) // '.h5'
// '_' // trim(to_str(p%id)) // '.h5'
!$omp critical (WriteParticleRestart)
! Create file
call pr % file_create(filename)
file_id = file_create(filename)
! Get information about source particle
src => source_bank(current_work)
! Write data to file
call pr % write_data(FILETYPE_PARTICLE_RESTART, 'filetype')
call pr % write_data(REVISION_PARTICLE_RESTART, 'revision')
call pr % write_data(current_batch, 'current_batch')
call pr % write_data(gen_per_batch, 'gen_per_batch')
call pr % write_data(current_gen, 'current_gen')
call pr % write_data(n_particles, 'n_particles')
call pr % write_data(run_mode, 'run_mode')
call pr % write_data(p % id, 'id')
call pr % write_data(src % wgt, 'weight')
call pr % write_data(src % E, 'energy')
call pr % write_data(src % xyz, 'xyz', length = 3)
call pr % write_data(src % uvw, 'uvw', length = 3)
call write_dataset(file_id, 'filetype', FILETYPE_PARTICLE_RESTART)
call write_dataset(file_id, 'revision', REVISION_PARTICLE_RESTART)
call write_dataset(file_id, 'current_batch', current_batch)
call write_dataset(file_id, 'gen_per_batch', gen_per_batch)
call write_dataset(file_id, 'current_gen', current_gen)
call write_dataset(file_id, 'n_particles', n_particles)
call write_dataset(file_id, 'run_mode', run_mode)
call write_dataset(file_id, 'id', p%id)
call write_dataset(file_id, 'weight', src%wgt)
call write_dataset(file_id, 'energy', src%E)
call write_dataset(file_id, 'xyz', src%xyz)
call write_dataset(file_id, 'uvw', src%uvw)
! Close file
call pr % file_close()
call file_close(file_id)
!$omp end critical (WriteParticleRestart)
end subroutine write_particle_restart

93
src/source.F90
View file

@ -6,17 +6,20 @@ module source
use geometry, only: find_cell
use geometry_header, only: BASE_UNIVERSE
use global
use hdf5_interface, only: file_create, file_open, file_close, read_dataset
use math, only: maxwell_spectrum, watt_spectrum
use output, only: write_message
use output_interface, only: BinaryOutput
use particle_header, only: Particle
use random_lcg, only: prn, set_particle_seed, prn_set_stream
use state_point, only: read_source_bank, write_source_bank
use string, only: to_str
#ifdef MPI
use message_passing
#endif
use hdf5, only: HID_T
implicit none
contains
@ -27,12 +30,12 @@ contains
subroutine initialize_source()
character(MAX_FILE_LEN) :: filename
integer(8) :: i ! loop index over bank sites
integer(8) :: id ! particle id
integer(4) :: itmp ! temporary integer
type(Bank), pointer :: src => null() ! source bank site
type(BinaryOutput) :: sp ! statepoint/source binary file
integer(HID_T) :: file_id
character(MAX_FILE_LEN) :: filename
type(Bank), pointer :: src ! source bank site
call write_message("Initializing source particles...", 6)
@ -44,10 +47,10 @@ contains
&// '...', 6)
! Open the binary file
call sp % file_open(path_source, 'r', serial = .false.)
file_id = file_open(path_source, 'r', parallel=.true.)
! Read the file type
call sp % read_data(itmp, "filetype")
call read_dataset(file_id, "filetype", itmp)
! Check to make sure this is a source file
if (itmp /= FILETYPE_SOURCE) then
@ -56,10 +59,10 @@ contains
end if
! Read in the source bank
call sp % read_source_bank()
call read_source_bank(file_id)
! Close file
call sp % file_close()
call file_close(file_id)
else
! Generation source sites from specified distribution in user input
@ -80,9 +83,9 @@ contains
if (write_initial_source) then
call write_message('Writing out initial source...', 1)
filename = trim(path_output) // 'initial_source.h5'
call sp % file_create(filename, serial = .false.)
call sp % write_source_bank()
call sp % file_close()
file_id = file_create(filename, parallel=.true.)
call write_source_bank(file_id)
call file_close(file_id)
end if
end subroutine initialize_source
@ -109,28 +112,28 @@ contains
integer, save :: num_resamples = 0 ! Number of resamples encountered
! Set weight to one by default
site % wgt = ONE
site%wgt = ONE
! Set the random number generator to the source stream.
call prn_set_stream(STREAM_SOURCE)
! Sample position
select case (external_source % type_space)
select case (external_source%type_space)
case (SRC_SPACE_BOX)
! Set particle defaults
call p % initialize()
call p%initialize()
! Repeat sampling source location until a good site has been found
found = .false.
do while (.not.found)
! Coordinates sampled uniformly over a box
p_min = external_source % params_space(1:3)
p_max = external_source % params_space(4:6)
p_min = external_source%params_space(1:3)
p_max = external_source%params_space(4:6)
r = (/ (prn(), i = 1,3) /)
site % xyz = p_min + r*(p_max - p_min)
site%xyz = p_min + r*(p_max - p_min)
! Fill p with needed data
p % coord(1) % xyz = site % xyz
p % coord(1) % uvw = [ ONE, ZERO, ZERO ]
p%coord(1)%xyz = site%xyz
p%coord(1)%uvw = [ ONE, ZERO, ZERO ]
! Now search to see if location exists in geometry
call find_cell(p, found)
@ -142,24 +145,24 @@ contains
end if
end if
end do
call p % clear()
call p%clear()
case (SRC_SPACE_FISSION)
! Repeat sampling source location until a good site has been found
found = .false.
do while (.not.found)
! Set particle defaults
call p % initialize()
call p%initialize()
! Coordinates sampled uniformly over a box
p_min = external_source % params_space(1:3)
p_max = external_source % params_space(4:6)
p_min = external_source%params_space(1:3)
p_max = external_source%params_space(4:6)
r = (/ (prn(), i = 1,3) /)
site % xyz = p_min + r*(p_max - p_min)
site%xyz = p_min + r*(p_max - p_min)
! Fill p with needed data
p % coord(1) % xyz = site % xyz
p % coord(1) % uvw = [ ONE, ZERO, ZERO ]
p%coord(1)%xyz = site%xyz
p%coord(1)%uvw = [ ONE, ZERO, ZERO ]
! Now search to see if location exists in geometry
call find_cell(p, found)
@ -171,66 +174,66 @@ contains
end if
cycle
end if
if (p % material == MATERIAL_VOID) then
if (p%material == MATERIAL_VOID) then
found = .false.
cycle
end if
if (.not. materials(p % material) % fissionable) found = .false.
if (.not. materials(p%material)%fissionable) found = .false.
end do
call p % clear()
call p%clear()
case (SRC_SPACE_POINT)
! Point source
site % xyz = external_source % params_space
site%xyz = external_source%params_space
end select
! Sample angle
select case (external_source % type_angle)
select case (external_source%type_angle)
case (SRC_ANGLE_ISOTROPIC)
! Sample isotropic distribution
phi = TWO*PI*prn()
mu = TWO*prn() - ONE
site % uvw(1) = mu
site % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
site % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
site%uvw(1) = mu
site%uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
site%uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
case (SRC_ANGLE_MONO)
! Monodirectional source
site % uvw = external_source % params_angle
site%uvw = external_source%params_angle
case default
call fatal_error("No angle distribution specified for external source!")
end select
! Sample energy distribution
select case (external_source % type_energy)
select case (external_source%type_energy)
case (SRC_ENERGY_MONO)
! Monoenergtic source
site % E = external_source % params_energy(1)
if (site % E >= 20) then
site%E = external_source%params_energy(1)
if (site%E >= 20) then
call fatal_error("Source energies above 20 MeV not allowed.")
end if
case (SRC_ENERGY_MAXWELL)
a = external_source % params_energy(1)
a = external_source%params_energy(1)
do
! Sample Maxwellian fission spectrum
site % E = maxwell_spectrum(a)
site%E = maxwell_spectrum(a)
! resample if energy is >= 20 MeV
if (site % E < 20) exit
if (site%E < 20) exit
end do
case (SRC_ENERGY_WATT)
a = external_source % params_energy(1)
b = external_source % params_energy(2)
a = external_source%params_energy(1)
b = external_source%params_energy(2)
do
! Sample Watt fission spectrum
site % E = watt_spectrum(a, b)
site%E = watt_spectrum(a, b)
! resample if energy is >= 20 MeV
if (site % E < 20) exit
if (site%E < 20) exit
end do
case default

828
src/state_point.F90
View file

File diff suppressed because it is too large Load diff

53
src/track_output.F90
View file

@ -6,19 +6,27 @@
module track_output
use global
use output_interface, only: BinaryOutput
use hdf5_interface
use particle_header, only: Particle
use string, only: to_str
implicit none
use hdf5
type, private :: TrackCoordinates
implicit none
private
type TrackCoordinates
real(8), allocatable :: coords(:,:)
end type TrackCoordinates
type(TrackCoordinates), private, allocatable :: tracks(:)
type(TrackCoordinates), allocatable :: tracks(:)
!$omp threadprivate(tracks)
public :: initialize_particle_track
public :: write_particle_track
public :: add_particle_track
public :: finalize_particle_track
contains
!===============================================================================
@ -43,19 +51,19 @@ contains
! Add another column to coords
i = size(tracks)
if (allocated(tracks(i) % coords)) then
n_tracks = size(tracks(i) % coords, 2)
if (allocated(tracks(i)%coords)) then
n_tracks = size(tracks(i)%coords, 2)
allocate(new_coords(3, n_tracks + 1))
new_coords(:, 1:n_tracks) = tracks(i) % coords
call move_alloc(FROM=new_coords, TO=tracks(i) % coords)
new_coords(:, 1:n_tracks) = tracks(i)%coords
call move_alloc(FROM=new_coords, TO=tracks(i)%coords)
else
n_tracks = 0
allocate(tracks(i) % coords(3, 1))
allocate(tracks(i)%coords(3, 1))
end if
! Write current coordinates into the newest column.
n_tracks = n_tracks + 1
tracks(i) % coords(:, n_tracks) = p % coord(1) % xyz
tracks(i)%coords(:, n_tracks) = p%coord(1)%xyz
end subroutine write_particle_track
!===============================================================================
@ -87,35 +95,32 @@ contains
subroutine finalize_particle_track(p)
type(Particle), intent(in) :: p
integer :: length(2)
character(MAX_FILE_LEN) :: fname
type(BinaryOutput) :: binout
integer :: i
integer, allocatable :: n_coords(:)
integer :: n_particle_tracks
integer(HID_T) :: file_id
character(MAX_FILE_LEN) :: fname
integer, allocatable :: n_coords(:)
fname = trim(path_output) // 'track_' // trim(to_str(current_batch)) &
// '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) &
// '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p%id)) &
// '.h5'
! Determine total number of particles and number of coordinates for each
n_particle_tracks = size(tracks)
allocate(n_coords(n_particle_tracks))
do i = 1, n_particle_tracks
n_coords(i) = size(tracks(i) % coords, 2)
n_coords(i) = size(tracks(i)%coords, 2)
end do
!$omp critical (FinalizeParticleTrack)
call binout % file_create(fname)
call binout % write_data(n_particle_tracks, 'n_particles')
call binout % write_data(n_coords, 'n_coords', length=n_particle_tracks)
file_id = file_create(fname)
call write_dataset(file_id, 'n_particles', n_particle_tracks)
call write_dataset(file_id, 'n_coords', n_coords)
do i = 1, n_particle_tracks
length(:) = [3, n_coords(i)]
call binout % write_data(tracks(i) % coords, 'coordinates_' // &
trim(to_str(i)), length=length)
call write_dataset(file_id, 'coordinates_' // trim(to_str(i)), &
tracks(i)%coords)
end do
call binout % file_close()
call file_close(file_id)
!$omp end critical (FinalizeParticleTrack)
deallocate(tracks)
end subroutine finalize_particle_track