From 4bfc53db159fbe7bbb3f01e161d0a10d012877fe Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 2 Apr 2013 05:16:38 -0700 Subject: [PATCH 01/21] added base file for particle restarts --- src/DEPENDENCIES | 4 +++ src/OBJECTS | 1 + src/particle_restart.F90 | 56 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 61 insertions(+) create mode 100644 src/particle_restart.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 61c31dd90..433f62f67 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -264,6 +264,10 @@ output.o: tally_header.o particle_header.o: constants.o +particle_restart.o: constants.o +particle_restart.o: global.o +particle_restart.o: hdf5_interface.o + physics.o: ace_header.o physics.o: constants.o physics.o: cross_section.o diff --git a/src/OBJECTS b/src/OBJECTS index d53f0a543..d9ec550f0 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -40,6 +40,7 @@ mesh_header.o \ mesh.o \ output.o \ particle_header.o \ +particle_restart.o \ physics.o \ plot.o \ plot_header.o \ diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 new file mode 100644 index 000000000..49cc5f3df --- /dev/null +++ b/src/particle_restart.F90 @@ -0,0 +1,56 @@ +module particle_restart + + use constants + use global + use string + +#ifdef HDF5 + use hdf5_interface +#endif + + private + public :: write_particle_restart + + integer(HID_T) :: hdf5_particle_file + +contains + +!=============================================================================== +! WRITE_PARTICLE_RESTART +!=============================================================================== + + subroutine write_particle_restart() + + character(MAX_FILE_LEN) :: filename + + ! set up file name + filename = 'particle_'//to_str(p % id)//'.h5' + + ! create hdf5 file + call h5fcreate_f(filename, H5F_ACC_TRUNC_F, hdf5_particle_file, hdf5_err) + + ! write data to file + + + ! close hdf5 file + call h5fclose_f(hdf5_particle_file, hdf5_err) + + end subroutine write_particle_restart + +!=============================================================================== +! READ_PARTICLE_RESTART +!=============================================================================== + + subroutine read_particle_restart + + end subroutine read_particle_restart + +!=============================================================================== +! RUN_PARTICLE_RESTART +!=============================================================================== + + subroutine run_particle_restart + + end subroutine run_particle_restart + +end module particle_restart From a95b48aa0eed70d3391904199f51801d9a9d90cf Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 2 Apr 2013 05:24:20 -0700 Subject: [PATCH 02/21] changed particle work counter local variable i to a global variable current_work --- src/eigenvalue.F90 | 4 ++-- src/global.F90 | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 8357c89f3..9c0678f7a 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -74,10 +74,10 @@ contains ! ==================================================================== ! LOOP OVER PARTICLES - PARTICLE_LOOP: do i = 1, work + PARTICLE_LOOP: do current_work = 1, work ! grab source particle from bank - call get_source_particle(i) + call get_source_particle(current_work) ! transport particle call transport() diff --git a/src/global.F90 b/src/global.F90 index 87e07de26..694241589 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -166,6 +166,7 @@ module global integer(8) :: bank_last ! index of last particle in bank integer(8) :: work ! number of particles per processor integer(8) :: maxwork ! maximum number of particles per processor + integer :: current_work ! current work ! Temporary k-effective values real(8), allocatable :: k_batch(:) ! batch estimates of k From 95b9d79db0470f6b13079eca3038c1c73d79d411 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 2 Apr 2013 06:07:31 -0700 Subject: [PATCH 03/21] changed current_work to a long integer --- src/global.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/global.F90 b/src/global.F90 index 694241589..fae7da8f1 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -166,7 +166,7 @@ module global integer(8) :: bank_last ! index of last particle in bank integer(8) :: work ! number of particles per processor integer(8) :: maxwork ! maximum number of particles per processor - integer :: current_work ! current work + integer(8) :: current_work ! current work ! Temporary k-effective values real(8), allocatable :: k_batch(:) ! batch estimates of k From 98a4105a5942a5bc316638001e1e2fbca0d40767 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 2 Apr 2013 06:32:11 -0700 Subject: [PATCH 04/21] particle restart file is now written from fatal error routine --- src/DEPENDENCIES | 2 +- src/error.F90 | 4 + src/particle_restart.F90 | 169 ++++++++++++++++++++++++++++++++++++++- 3 files changed, 170 insertions(+), 5 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 433f62f67..3e4d77b8a 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -123,6 +123,7 @@ energy_grid.o: list_header.o energy_grid.o: output.o error.o: global.o +error.o: particle_restart.o finalize.o: cmfd_output.o finalize.o: global.o @@ -266,7 +267,6 @@ particle_header.o: constants.o particle_restart.o: constants.o particle_restart.o: global.o -particle_restart.o: hdf5_interface.o physics.o: ace_header.o physics.o: constants.o diff --git a/src/error.F90 b/src/error.F90 index ec8b9f2fa..9d25ae655 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -3,6 +3,7 @@ module error use, intrinsic :: ISO_FORTRAN_ENV use global + use particle_restart, only: write_particle_restart #ifdef MPI use mpi @@ -145,6 +146,9 @@ contains write(ERROR_UNIT,*) end if + ! Write particle restart + call write_particle_restart() + ! Release memory from all allocatable arrays call free_memory() diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 49cc5f3df..da8b71b96 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -1,11 +1,13 @@ module particle_restart + use bank_header, only: Bank use constants use global - use string #ifdef HDF5 - use hdf5_interface + use hdf5 + use h5lt + use, intrinsic :: ISO_C_BINDING #endif private @@ -22,15 +24,31 @@ contains subroutine write_particle_restart() character(MAX_FILE_LEN) :: filename + integer(HSIZE_T) :: dims1(1) + type(Bank), pointer :: src => null() ! set up file name - filename = 'particle_'//to_str(p % id)//'.h5' + filename = 'particle_'//trim(int4_to_str(rank))//'.h5' ! create hdf5 file call h5fcreate_f(filename, H5F_ACC_TRUNC_F, hdf5_particle_file, hdf5_err) - ! write data to file + ! get information about source particle + src => source_bank(current_work) + ! write data to file + call hdf5_write_integer(hdf5_particle_file, 'current_batch', current_batch) + call hdf5_write_integer(hdf5_particle_file, 'gen_per_batch', gen_per_batch) + call hdf5_write_integer(hdf5_particle_file, 'current_gen', current_gen) + call hdf5_write_long(hdf5_particle_file, 'n_particles', n_particles) + call hdf5_write_long(hdf5_particle_file, 'id', p % id) + call hdf5_write_double(hdf5_particle_file, 'weight', src % wgt) + call hdf5_write_double(hdf5_particle_file, 'energy', src % E) + dims1 = (/3/) + call h5ltmake_dataset_double_f(hdf5_particle_file, 'xyz', 1, dims1, & + src % xyz, hdf5_err) + call h5ltmake_dataset_double_f(hdf5_particle_file, 'uvw', 1, dims1, & + src % uvw, hdf5_err) ! close hdf5 file call h5fclose_f(hdf5_particle_file, hdf5_err) @@ -53,4 +71,147 @@ contains end subroutine run_particle_restart +!=============================================================================== +! HDF5_WRITE_INTEGER +!=============================================================================== + + subroutine hdf5_write_integer(group, name, buffer) + + integer(HID_T), intent(in) :: group + character(*), intent(in) :: name + integer, intent(in) :: buffer + + integer :: rank = 1 + integer(HSIZE_T) :: dims(1) = (/1/) + + call h5ltmake_dataset_int_f(group, name, rank, dims, & + (/ buffer /), hdf5_err) + + end subroutine hdf5_write_integer + +!=============================================================================== +! HDF5_WRITE_LONG +!=============================================================================== + + subroutine hdf5_write_long(group, name, buffer) + + integer(HID_T), intent(in) :: group + character(*), intent(in) :: name + integer(8), target, intent(in) :: buffer + + integer :: rank = 1 + integer(HSIZE_T) :: dims(1) = (/1/) + integer(HID_T) :: dspace + integer(HID_T) :: dset + type(c_ptr) :: f_ptr + + ! Create dataspace and dataset + call h5screate_simple_f(rank, dims, dspace, hdf5_err) + call h5dcreate_f(group, name, hdf5_integer8_t, dspace, dset, hdf5_err) + + ! Write eight-byte integer + f_ptr = c_loc(buffer) + call h5dwrite_f(dset, hdf5_integer8_t, f_ptr, hdf5_err) + + ! Close dataspace and dataset for long integer + call h5dclose_f(dset, hdf5_err) + call h5sclose_f(dspace, hdf5_err) + + end subroutine hdf5_write_long + +!=============================================================================== +! HDF5_WRITE_DOUBLE +!=============================================================================== + + subroutine hdf5_write_double(group, name, buffer) + + integer(HID_T), intent(in) :: group + character(*), intent(in) :: name + real(8), intent(in) :: buffer + + integer :: rank = 1 + integer(HSIZE_T) :: dims(1) = (/1/) + + call h5ltmake_dataset_double_f(group, name, rank, dims, & + (/ buffer /), hdf5_err) + + end subroutine hdf5_write_double + +!=============================================================================== +! HDF5_READ_INTEGER +!=============================================================================== + + subroutine hdf5_read_integer(group, name, buffer) + + integer(HID_T), intent(in) :: group + character(*), intent(in) :: name + integer, intent(inout) :: buffer + + integer :: buffer_copy(1) + integer(HSIZE_T) :: dims(1) = (/1/) + + call h5ltread_dataset_int_f(group, name, buffer_copy, dims, hdf5_err) + buffer = buffer_copy(1) + + end subroutine hdf5_read_integer + +!=============================================================================== +! HDF5_READ_LONG +!=============================================================================== + + subroutine hdf5_read_long(group, name, buffer) + + integer(HID_T), intent(in) :: group + character(*), intent(in) :: name + integer(8), target, intent(out) :: buffer + + integer(HID_T) :: dset + type(c_ptr) :: f_ptr + + ! Open dataset + call h5dopen_f(group, name, dset, hdf5_err) + + ! Get pointer to buffer + f_ptr = c_loc(buffer) + + ! Read data from dataset + call h5dread_f(dset, hdf5_integer8_t, f_ptr, hdf5_err) + + ! Close dataset + call h5dclose_f(dset, hdf5_err) + + end subroutine hdf5_read_long + +!=============================================================================== +! HDF5_READ_DOUBLE +!=============================================================================== + + subroutine hdf5_read_double(group, name, buffer) + + integer(HID_T), intent(in) :: group + character(*), intent(in) :: name + real(8), intent(out) :: buffer + + real(8) :: buffer_copy(1) + integer(HSIZE_T) :: dims(1) = (/1/) + + call h5ltread_dataset_double_f(group, name, buffer_copy, dims, hdf5_err) + buffer = buffer_copy(1) + + end subroutine hdf5_read_double + +!=============================================================================== +! INT4_TO_STR converts an integer(4) to a string. +!=============================================================================== + + function int4_to_str(num) result(str) + + integer, intent(in) :: num + character(11) :: str + + write (str, '(I11)') num + str = adjustl(str) + + end function int4_to_str + end module particle_restart From 42c1159aa5e83811628c0ddac43048d2634ae305 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 2 Apr 2013 07:51:54 -0700 Subject: [PATCH 05/21] added new run mode for particle restart, restart files are read in --- src/DEPENDENCIES | 3 ++ src/constants.F90 | 3 +- src/finalize.F90 | 5 ++- src/global.F90 | 12 ++++-- src/initialize.F90 | 8 ++++ src/main.F90 | 13 +++--- src/particle_restart.F90 | 86 +++++++++++++++++++++++++++++++++++++--- 7 files changed, 113 insertions(+), 17 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 3e4d77b8a..11c2b1370 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -265,8 +265,11 @@ output.o: tally_header.o particle_header.o: constants.o +particle_restart.o: bank_header.o particle_restart.o: constants.o +particle_restart.o: geometry_header.o particle_restart.o: global.o +particle_restart.o: particle_header.o physics.o: ace_header.o physics.o: constants.o diff --git a/src/constants.F90 b/src/constants.F90 index 1a3d330ff..156f0046a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -346,7 +346,8 @@ module constants MODE_FIXEDSOURCE = 1, & ! Fixed source mode MODE_EIGENVALUE = 2, & ! K eigenvalue mode MODE_PLOTTING = 3, & ! Plotting mode - MODE_TALLIES = 4 ! Tally results mode + MODE_TALLIES = 4, & ! Tally results mode + MODE_PARTICLE = 5 ! Particle restart mode ! Unit numbers integer, parameter :: UNIT_SUMMARY = 11 ! unit # for writing summary file diff --git a/src/finalize.F90 b/src/finalize.F90 index 1bfdfe71f..d601c2e59 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -29,7 +29,7 @@ contains ! Start finalization timer call time_finalize % start() - if (run_mode /= MODE_PLOTTING) then + if (run_mode /= MODE_PLOTTING .and. run_mode /= MODE_PARTICLE) then ! Calculate statistics for tallies and write to tallies.out if (master) call tally_statistics() if (output_tallies) then @@ -46,7 +46,8 @@ contains call time_finalize % stop() call time_total % stop() if (master .and. (run_mode /= MODE_PLOTTING .and. & - run_mode /= MODE_TALLIES)) then + run_mode /= MODE_TALLIES .and. & + run_mode /= MODE_PARTICLE)) then call print_runtime() call print_results() end if diff --git a/src/global.F90 b/src/global.F90 index fae7da8f1..12c88fc58 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -254,10 +254,11 @@ module global logical :: restart_run = .false. integer :: restart_batch - character(MAX_FILE_LEN) :: path_input ! Path to input file - character(MAX_FILE_LEN) :: path_cross_sections ! Path to cross_sections.xml - character(MAX_FILE_LEN) :: path_source = '' ! Path to binary source - character(MAX_FILE_LEN) :: path_state_point ! Path to binary state point + character(MAX_FILE_LEN) :: path_input ! Path to input file + character(MAX_FILE_LEN) :: path_cross_sections ! Path to cross_sections.xml + character(MAX_FILE_LEN) :: path_source = '' ! Path to binary source + character(MAX_FILE_LEN) :: path_state_point ! Path to binary state point + character(MAX_FILE_LEN) :: path_particle_restart ! Path to particle restart ! Message used in message/warning/fatal_error character(MAX_LINE_LEN) :: message @@ -275,6 +276,9 @@ module global integer :: trace_gen integer(8) :: trace_particle + ! Particle restart run + logical :: particle_restart_run = .false. + ! ============================================================================ ! CMFD VARIABLES diff --git a/src/initialize.F90 b/src/initialize.F90 index 944f84274..33fb953d5 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -145,6 +145,9 @@ contains end if end if + ! check for particle restart run + if (particle_restart_run) run_mode = MODE_PARTICLE + ! Stop initialization timer call time_initialize % stop() @@ -303,6 +306,11 @@ contains case ('-eps_tol', '-ksp_gmres_restart') ! Handle options that would be based to PETSC i = i + 1 + case ('-partcle','--particle') + ! Read in path for particle restart + i = i + 1 + path_particle_restart = argv(i) + particle_restart_run = .true. case default message = "Unknown command line option: " // argv(i) call fatal_error() diff --git a/src/main.F90 b/src/main.F90 index 584db7f1e..41a7481da 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -1,12 +1,13 @@ program main use constants - use eigenvalue, only: run_eigenvalue - use finalize, only: finalize_run - use fixed_source, only: run_fixedsource + use eigenvalue, only: run_eigenvalue + use finalize, only: finalize_run + use fixed_source, only: run_fixedsource use global - use initialize, only: initialize_run - use plot, only: run_plot + use initialize, only: initialize_run + use particle_restart, only: run_particle_restart + use plot, only: run_plot implicit none @@ -24,6 +25,8 @@ program main case (MODE_TALLIES) ! For tallies-only mode, we just skip straight to finalize_run to write out ! the tally results + case (MODE_PARTICLE) + call run_particle_restart() end select ! finalize run diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index da8b71b96..2e08c9c93 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -1,17 +1,19 @@ module particle_restart - use bank_header, only: Bank + use bank_header, only: Bank use constants + use geometry_header, only: BASE_UNIVERSE use global + use particle_header, only: deallocate_coord #ifdef HDF5 use hdf5 use h5lt - use, intrinsic :: ISO_C_BINDING #endif + implicit none private - public :: write_particle_restart + public :: write_particle_restart, run_particle_restart integer(HID_T) :: hdf5_particle_file @@ -59,7 +61,35 @@ contains ! READ_PARTICLE_RESTART !=============================================================================== - subroutine read_particle_restart + subroutine read_particle_restart() + + integer(HSIZE_T) :: dims1(1) + + ! open hdf5 file + call h5fopen_f(path_particle_restart, H5F_ACC_RDONLY_F, hdf5_particle_file,& + hdf5_err) + + ! read data from file + call hdf5_read_integer(hdf5_particle_file, 'current_batch', current_batch) + call hdf5_read_integer(hdf5_particle_file, 'gen_per_batch', gen_per_batch) + call hdf5_read_integer(hdf5_particle_file, 'current_gen', current_gen) + call hdf5_read_long(hdf5_particle_file, 'n_particles', n_particles) + call hdf5_read_long(hdf5_particle_file, 'id', p % id) + call hdf5_read_double(hdf5_particle_file, 'weight', p % wgt) + call hdf5_read_double(hdf5_particle_file, 'energy', p % E) + dims1 = (/3/) + call h5ltread_dataset_double_f(hdf5_particle_file, 'xyz', p % coord % xyz, & + dims1, hdf5_err) + call h5ltread_dataset_double_f(hdf5_particle_file, 'uvw', p % coord % uvw, & + dims1, hdf5_err) + + ! set particle last attributes + p % last_wgt = p % wgt + p % last_xyz = p % coord % xyz + p % last_E = p % E + + ! close hdf5 file + call h5fclose_f(hdf5_particle_file, hdf5_err) end subroutine read_particle_restart @@ -67,10 +97,56 @@ contains ! RUN_PARTICLE_RESTART !=============================================================================== - subroutine run_particle_restart + subroutine run_particle_restart() + + ! initialize the particle to be tracked + call initialize_particle() + + ! read in the restart information + call read_particle_restart() + + ! set all tallies to 0 for now (just tracking errors) + n_tallies = 0 + + ! compute seed end subroutine run_particle_restart +!=============================================================================== +! INITIALIZE_PARTICLE +!=============================================================================== + + subroutine initialize_particle() + + ! Allocate particle + allocate(p) + + ! Set particle to neutron that's alive + p % type = NEUTRON + p % alive = .true. + + ! clear attributes + p % surface = NONE + p % cell_born = NONE + p % material = NONE + p % last_material = NONE + p % wgt = ONE + p % last_wgt = ONE + p % absorb_wgt = ZERO + p % n_bank = 0 + p % wgt_bank = ZERO + p % n_collision = 0 + + ! remove any original coordinates + call deallocate_coord(p % coord0) + + ! Set up base level coordinates + allocate(p % coord0) + p % coord0 % universe = BASE_UNIVERSE + p % coord => p % coord0 + + end subroutine initialize_particle + !=============================================================================== ! HDF5_WRITE_INTEGER !=============================================================================== From 221231f82cc8f9f348a93bccba91bb29cba13db7 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 2 Apr 2013 08:26:28 -0700 Subject: [PATCH 06/21] added HDF5 compiler directives and added output to user --- src/particle_restart.F90 | 82 +++++++++++++++++++++++++++++++++++----- 1 file changed, 72 insertions(+), 10 deletions(-) diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 2e08c9c93..bc59d6677 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -1,5 +1,7 @@ module particle_restart + use, intrinsic :: ISO_FORTRAN_ENV + use bank_header, only: Bank use constants use geometry_header, only: BASE_UNIVERSE @@ -15,7 +17,13 @@ module particle_restart private public :: write_particle_restart, run_particle_restart +#ifdef HDF5 integer(HID_T) :: hdf5_particle_file +#endif HDF5 + + ! Short names for output and error units + integer :: ou = OUTPUT_UNIT + integer :: eu = ERROR_UNIT contains @@ -24,7 +32,7 @@ contains !=============================================================================== subroutine write_particle_restart() - +#ifdef HDF5 character(MAX_FILE_LEN) :: filename integer(HSIZE_T) :: dims1(1) type(Bank), pointer :: src => null() @@ -54,17 +62,22 @@ contains ! close hdf5 file call h5fclose_f(hdf5_particle_file, hdf5_err) - +#endif end subroutine write_particle_restart !=============================================================================== ! READ_PARTICLE_RESTART !=============================================================================== - subroutine read_particle_restart() - + subroutine read_hdf5_particle_restart() +#ifdef HDF5 integer(HSIZE_T) :: dims1(1) + ! write meessage + message = "Loading particle restart file " // trim(path_particle_restart) & + // "..." + call write_message(1) + ! open hdf5 file call h5fopen_f(path_particle_restart, H5F_ACC_RDONLY_F, hdf5_particle_file,& hdf5_err) @@ -90,26 +103,26 @@ contains ! close hdf5 file call h5fclose_f(hdf5_particle_file, hdf5_err) - - end subroutine read_particle_restart +#endif + end subroutine read_hdf5_particle_restart !=============================================================================== ! RUN_PARTICLE_RESTART !=============================================================================== subroutine run_particle_restart() - +#ifdef HDF5 ! initialize the particle to be tracked call initialize_particle() ! read in the restart information - call read_particle_restart() + call read_hdf5_particle_restart() ! set all tallies to 0 for now (just tracking errors) n_tallies = 0 ! compute seed - +#endif HDF5 end subroutine run_particle_restart !=============================================================================== @@ -150,7 +163,7 @@ contains !=============================================================================== ! HDF5_WRITE_INTEGER !=============================================================================== - +#ifdef HDF5 subroutine hdf5_write_integer(group, name, buffer) integer(HID_T), intent(in) :: group @@ -289,5 +302,54 @@ contains str = adjustl(str) end function int4_to_str +#endif + +!=============================================================================== +! WRITE_MESSAGE displays an informational message to the log file and the +! standard output stream. +!=============================================================================== + + subroutine write_message(level) + + integer, optional :: level ! verbosity level + + integer :: i_start ! starting position + integer :: i_end ! ending position + integer :: line_wrap ! length of line + integer :: length ! length of message + + ! Set length of line + line_wrap = 80 + + ! Only allow master to print to screen + if (.not. master .and. present(level)) return + + if (.not. present(level) .or. level <= verbosity) then + ! Determine length of message + length = len_trim(message) + + i_start = 0 + do + if (length - i_start < line_wrap - 1) then + ! Remainder of message will fit on line + write(ou, fmt='(1X,A)') message(i_start+1:length) + exit + + else + ! Determine last space in current line + i_end = i_start + index(message(i_start+1:i_start+line_wrap), & + ' ', BACK=.true.) + + ! Write up to last space + write(ou, fmt='(1X,A)') message(i_start+1:i_end-1) + + ! Advance starting position + i_start = i_end + if (i_start > length) exit + end if + end do + end if + + end subroutine write_message end module particle_restart From 6134642f35142d8fc1bc1824d35c95304b8b5691 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 2 Apr 2013 09:45:03 -0700 Subject: [PATCH 07/21] split particle restart logic into two files, call to write particle restart is in physics and geometry now decoupled from error --- src/DEPENDENCIES | 14 +- src/OBJECTS | 1 + src/error.F90 | 4 - src/geometry.F90 | 19 ++- src/particle_restart.F90 | 291 +++------------------------------ src/particle_restart_write.F90 | 61 +++++++ src/physics.F90 | 62 +++++-- 7 files changed, 154 insertions(+), 298 deletions(-) create mode 100644 src/particle_restart_write.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 11c2b1370..5e2a96a09 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -123,7 +123,6 @@ energy_grid.o: list_header.o energy_grid.o: output.o error.o: global.o -error.o: particle_restart.o finalize.o: cmfd_output.o finalize.o: global.o @@ -155,6 +154,7 @@ geometry.o: geometry_header.o geometry.o: global.o geometry.o: output.o geometry.o: particle_header.o +geometry.o: particle_restart_write.o geometry.o: string.o geometry.o: tally.o @@ -240,6 +240,7 @@ main.o: finalize.o main.o: fixed_source.o main.o: global.o main.o: initialize.o +main.o: particle_restart.o main.o: plot.o math.o: constants.o @@ -269,7 +270,17 @@ particle_restart.o: bank_header.o particle_restart.o: constants.o particle_restart.o: geometry_header.o particle_restart.o: global.o +particle_restart.o: hdf5_interface.o particle_restart.o: particle_header.o +particle_restart.o: output.o +particle_restart.o: physics.o +particle_restart.o: random_lcg.o +particle_restart.o: source.o + +particle_restart_write.o: bank_header.o +particle_restart_write.o: global.o +particle_restart_write.o: hdf5_interface.o +particle_restart_write.o: string.o physics.o: ace_header.o physics.o: constants.o @@ -285,6 +296,7 @@ physics.o: material_header.o physics.o: mesh.o physics.o: output.o physics.o: particle_header.o +physics.o: particle_restart_write.o physics.o: random_lcg.o physics.o: search.o physics.o: string.o diff --git a/src/OBJECTS b/src/OBJECTS index d9ec550f0..c9612e6af 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -41,6 +41,7 @@ mesh.o \ output.o \ particle_header.o \ particle_restart.o \ +particle_restart_write.o \ physics.o \ plot.o \ plot_header.o \ diff --git a/src/error.F90 b/src/error.F90 index 9d25ae655..ec8b9f2fa 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -3,7 +3,6 @@ module error use, intrinsic :: ISO_FORTRAN_ENV use global - use particle_restart, only: write_particle_restart #ifdef MPI use mpi @@ -146,9 +145,6 @@ contains write(ERROR_UNIT,*) end if - ! Write particle restart - call write_particle_restart() - ! Release memory from all allocatable arrays call free_memory() diff --git a/src/geometry.F90 b/src/geometry.F90 index 64912cb95..3ab85ba72 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -1,13 +1,14 @@ module geometry use constants - use error, only: fatal_error - use geometry_header, only: Cell, Surface, Universe, Lattice + use error, only: fatal_error + use geometry_header, only: Cell, Surface, Universe, Lattice use global - use output, only: write_message - use particle_header, only: LocalCoord, deallocate_coord - use string, only: to_str - use tally, only: score_surface_current + use output, only: write_message + use particle_header, only: LocalCoord, deallocate_coord + use particle_restart_write, only: write_particle_restart + use string, only: to_str + use tally, only: score_surface_current implicit none @@ -303,6 +304,7 @@ contains ! Do not handle reflective boundary conditions on lower universes if (.not. associated(p % coord, p % coord0)) then + call write_particle_restart() message = "Cannot reflect particle " // trim(to_str(p % id)) // & " off surface in a lower universe." call fatal_error() @@ -436,6 +438,7 @@ contains w = w + 2*dot_prod*R*z case default + call write_particle_restart() message = "Reflection not supported for surface " // & trim(to_str(surf % id)) call fatal_error() @@ -456,6 +459,7 @@ contains call deallocate_coord(p % coord0 % next) call find_cell(found) if (.not. found) then + call write_particle_restart() message = "Couldn't find particle after reflecting from surface." call fatal_error() end if @@ -515,6 +519,7 @@ contains ! undefined region in the geometry. if (.not. found) then + call write_particle_restart() message = "After particle " // trim(to_str(p % id)) // " crossed surface " & // trim(to_str(surfaces(abs(p%surface)) % id)) // " it could not be & &located in any cell and it did not leak." @@ -608,6 +613,7 @@ contains ! Search for particle call find_cell(found) if (.not. found) then + call write_particle_restart() message = "Could not locate particle " // trim(to_str(p % id)) // & " after crossing a lattice boundary." call fatal_error() @@ -630,6 +636,7 @@ contains ! Search for particle call find_cell(found) if (.not. found) then + call write_particle_restart() message = "Could not locate particle " // trim(to_str(p % id)) // & " after crossing a lattice boundary." call fatal_error() diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index bc59d6677..ac2730b8f 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -7,15 +7,18 @@ module particle_restart use geometry_header, only: BASE_UNIVERSE use global use particle_header, only: deallocate_coord + use output, only: write_message + use physics, only: transport + use random_lcg, only: set_particle_seed + use source, only: initialize_particle #ifdef HDF5 - use hdf5 - use h5lt + use hdf5_interface #endif implicit none private - public :: write_particle_restart, run_particle_restart + public :: run_particle_restart #ifdef HDF5 integer(HID_T) :: hdf5_particle_file @@ -27,44 +30,6 @@ module particle_restart contains -!=============================================================================== -! WRITE_PARTICLE_RESTART -!=============================================================================== - - subroutine write_particle_restart() -#ifdef HDF5 - character(MAX_FILE_LEN) :: filename - integer(HSIZE_T) :: dims1(1) - type(Bank), pointer :: src => null() - - ! set up file name - filename = 'particle_'//trim(int4_to_str(rank))//'.h5' - - ! create hdf5 file - call h5fcreate_f(filename, H5F_ACC_TRUNC_F, hdf5_particle_file, hdf5_err) - - ! get information about source particle - src => source_bank(current_work) - - ! write data to file - call hdf5_write_integer(hdf5_particle_file, 'current_batch', current_batch) - call hdf5_write_integer(hdf5_particle_file, 'gen_per_batch', gen_per_batch) - call hdf5_write_integer(hdf5_particle_file, 'current_gen', current_gen) - call hdf5_write_long(hdf5_particle_file, 'n_particles', n_particles) - call hdf5_write_long(hdf5_particle_file, 'id', p % id) - call hdf5_write_double(hdf5_particle_file, 'weight', src % wgt) - call hdf5_write_double(hdf5_particle_file, 'energy', src % E) - dims1 = (/3/) - call h5ltmake_dataset_double_f(hdf5_particle_file, 'xyz', 1, dims1, & - src % xyz, hdf5_err) - call h5ltmake_dataset_double_f(hdf5_particle_file, 'uvw', 1, dims1, & - src % uvw, hdf5_err) - - ! close hdf5 file - call h5fclose_f(hdf5_particle_file, hdf5_err) -#endif - end subroutine write_particle_restart - !=============================================================================== ! READ_PARTICLE_RESTART !=============================================================================== @@ -112,7 +77,11 @@ contains subroutine run_particle_restart() #ifdef HDF5 + + integer(8) :: particle_seed + ! initialize the particle to be tracked + allocate(p) call initialize_particle() ! read in the restart information @@ -121,235 +90,19 @@ contains ! set all tallies to 0 for now (just tracking errors) n_tallies = 0 - ! compute seed + ! compute random number seed + particle_seed = ((current_batch - 1)*gen_per_batch + & + current_gen - 1)*n_particles + p % id + call set_particle_seed(particle_seed) + + ! transport neutron + call transport() + print *, 'WEIGHT:', p % wgt + print *, 'ENERGY:', p % E + print *, 'LOCATION:', p % coord % xyz + print *, 'ANGLE:', p % coord % uvw + #endif HDF5 end subroutine run_particle_restart -!=============================================================================== -! INITIALIZE_PARTICLE -!=============================================================================== - - subroutine initialize_particle() - - ! Allocate particle - allocate(p) - - ! Set particle to neutron that's alive - p % type = NEUTRON - p % alive = .true. - - ! clear attributes - p % surface = NONE - p % cell_born = NONE - p % material = NONE - p % last_material = NONE - p % wgt = ONE - p % last_wgt = ONE - p % absorb_wgt = ZERO - p % n_bank = 0 - p % wgt_bank = ZERO - p % n_collision = 0 - - ! remove any original coordinates - call deallocate_coord(p % coord0) - - ! Set up base level coordinates - allocate(p % coord0) - p % coord0 % universe = BASE_UNIVERSE - p % coord => p % coord0 - - end subroutine initialize_particle - -!=============================================================================== -! HDF5_WRITE_INTEGER -!=============================================================================== -#ifdef HDF5 - subroutine hdf5_write_integer(group, name, buffer) - - integer(HID_T), intent(in) :: group - character(*), intent(in) :: name - integer, intent(in) :: buffer - - integer :: rank = 1 - integer(HSIZE_T) :: dims(1) = (/1/) - - call h5ltmake_dataset_int_f(group, name, rank, dims, & - (/ buffer /), hdf5_err) - - end subroutine hdf5_write_integer - -!=============================================================================== -! HDF5_WRITE_LONG -!=============================================================================== - - subroutine hdf5_write_long(group, name, buffer) - - integer(HID_T), intent(in) :: group - character(*), intent(in) :: name - integer(8), target, intent(in) :: buffer - - integer :: rank = 1 - integer(HSIZE_T) :: dims(1) = (/1/) - integer(HID_T) :: dspace - integer(HID_T) :: dset - type(c_ptr) :: f_ptr - - ! Create dataspace and dataset - call h5screate_simple_f(rank, dims, dspace, hdf5_err) - call h5dcreate_f(group, name, hdf5_integer8_t, dspace, dset, hdf5_err) - - ! Write eight-byte integer - f_ptr = c_loc(buffer) - call h5dwrite_f(dset, hdf5_integer8_t, f_ptr, hdf5_err) - - ! Close dataspace and dataset for long integer - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - - end subroutine hdf5_write_long - -!=============================================================================== -! HDF5_WRITE_DOUBLE -!=============================================================================== - - subroutine hdf5_write_double(group, name, buffer) - - integer(HID_T), intent(in) :: group - character(*), intent(in) :: name - real(8), intent(in) :: buffer - - integer :: rank = 1 - integer(HSIZE_T) :: dims(1) = (/1/) - - call h5ltmake_dataset_double_f(group, name, rank, dims, & - (/ buffer /), hdf5_err) - - end subroutine hdf5_write_double - -!=============================================================================== -! HDF5_READ_INTEGER -!=============================================================================== - - subroutine hdf5_read_integer(group, name, buffer) - - integer(HID_T), intent(in) :: group - character(*), intent(in) :: name - integer, intent(inout) :: buffer - - integer :: buffer_copy(1) - integer(HSIZE_T) :: dims(1) = (/1/) - - call h5ltread_dataset_int_f(group, name, buffer_copy, dims, hdf5_err) - buffer = buffer_copy(1) - - end subroutine hdf5_read_integer - -!=============================================================================== -! HDF5_READ_LONG -!=============================================================================== - - subroutine hdf5_read_long(group, name, buffer) - - integer(HID_T), intent(in) :: group - character(*), intent(in) :: name - integer(8), target, intent(out) :: buffer - - integer(HID_T) :: dset - type(c_ptr) :: f_ptr - - ! Open dataset - call h5dopen_f(group, name, dset, hdf5_err) - - ! Get pointer to buffer - f_ptr = c_loc(buffer) - - ! Read data from dataset - call h5dread_f(dset, hdf5_integer8_t, f_ptr, hdf5_err) - - ! Close dataset - call h5dclose_f(dset, hdf5_err) - - end subroutine hdf5_read_long - -!=============================================================================== -! HDF5_READ_DOUBLE -!=============================================================================== - - subroutine hdf5_read_double(group, name, buffer) - - integer(HID_T), intent(in) :: group - character(*), intent(in) :: name - real(8), intent(out) :: buffer - - real(8) :: buffer_copy(1) - integer(HSIZE_T) :: dims(1) = (/1/) - - call h5ltread_dataset_double_f(group, name, buffer_copy, dims, hdf5_err) - buffer = buffer_copy(1) - - end subroutine hdf5_read_double - -!=============================================================================== -! INT4_TO_STR converts an integer(4) to a string. -!=============================================================================== - - function int4_to_str(num) result(str) - - integer, intent(in) :: num - character(11) :: str - - write (str, '(I11)') num - str = adjustl(str) - - end function int4_to_str -#endif - -!=============================================================================== -! WRITE_MESSAGE displays an informational message to the log file and the -! standard output stream. -!=============================================================================== - - subroutine write_message(level) - - integer, optional :: level ! verbosity level - - integer :: i_start ! starting position - integer :: i_end ! ending position - integer :: line_wrap ! length of line - integer :: length ! length of message - - ! Set length of line - line_wrap = 80 - - ! Only allow master to print to screen - if (.not. master .and. present(level)) return - - if (.not. present(level) .or. level <= verbosity) then - ! Determine length of message - length = len_trim(message) - - i_start = 0 - do - if (length - i_start < line_wrap - 1) then - ! Remainder of message will fit on line - write(ou, fmt='(1X,A)') message(i_start+1:length) - exit - - else - ! Determine last space in current line - i_end = i_start + index(message(i_start+1:i_start+line_wrap), & - ' ', BACK=.true.) - - ! Write up to last space - write(ou, fmt='(1X,A)') message(i_start+1:i_end-1) - - ! Advance starting position - i_start = i_end - if (i_start > length) exit - end if - end do - end if - - end subroutine write_message - end module particle_restart diff --git a/src/particle_restart_write.F90 b/src/particle_restart_write.F90 new file mode 100644 index 000000000..d9388ca76 --- /dev/null +++ b/src/particle_restart_write.F90 @@ -0,0 +1,61 @@ +module particle_restart_write + + use, intrinsic :: ISO_FORTRAN_ENV + + use bank_header, only: Bank + use global + use string, only: to_str + +#ifdef HDF5 + use hdf5_interface +#endif + + implicit none + private + public :: write_particle_restart + +#ifdef HDF5 + integer(HID_T) :: hdf5_particle_file +#endif HDF5 + +contains + +!=============================================================================== +! WRITE_PARTICLE_RESTART +!=============================================================================== + + subroutine write_particle_restart() +#ifdef HDF5 + character(MAX_FILE_LEN) :: filename + integer(HSIZE_T) :: dims1(1) + type(Bank), pointer :: src => null() + + ! set up file name + filename = 'particle_'//trim(to_str(rank))//'.h5' + + ! create hdf5 file + call h5fcreate_f(filename, H5F_ACC_TRUNC_F, hdf5_particle_file, hdf5_err) + + ! get information about source particle + src => source_bank(current_work) + + ! write data to file + call hdf5_write_integer(hdf5_particle_file, 'current_batch', current_batch) + call hdf5_write_integer(hdf5_particle_file, 'gen_per_batch', gen_per_batch) + call hdf5_write_integer(hdf5_particle_file, 'current_gen', current_gen) + call hdf5_write_long(hdf5_particle_file, 'n_particles', n_particles) + call hdf5_write_long(hdf5_particle_file, 'id', p % id) + call hdf5_write_double(hdf5_particle_file, 'weight', src % wgt) + call hdf5_write_double(hdf5_particle_file, 'energy', src % E) + dims1 = (/3/) + call h5ltmake_dataset_double_f(hdf5_particle_file, 'xyz', 1, dims1, & + src % xyz, hdf5_err) + call h5ltmake_dataset_double_f(hdf5_particle_file, 'uvw', 1, dims1, & + src % uvw, hdf5_err) + + ! close hdf5 file + call h5fclose_f(hdf5_particle_file, hdf5_err) +#endif + end subroutine write_particle_restart + +end module particle_restart_write diff --git a/src/physics.F90 b/src/physics.F90 index b252597ee..76c4141cd 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,25 +1,27 @@ module physics - use ace_header, only: Nuclide, Reaction, DistEnergy + use ace_header, only: Nuclide, Reaction, DistEnergy use constants - use cross_section, only: calculate_xs - use endf, only: reaction_name - use error, only: fatal_error, warning - use fission, only: nu_total, nu_delayed - use geometry, only: find_cell, distance_to_boundary, cross_surface, & - cross_lattice - use geometry_header, only: Universe, BASE_UNIVERSE + use cross_section, only: calculate_xs + use endf, only: reaction_name + use error, only: fatal_error, warning + use fission, only: nu_total, nu_delayed + use geometry, only: find_cell, distance_to_boundary, & + cross_surface, cross_lattice + use geometry_header, only: Universe, BASE_UNIVERSE use global - use interpolation, only: interpolate_tab1 - use material_header, only: Material - use mesh, only: get_mesh_indices - use output, only: write_message - use particle_header, only: LocalCoord - use random_lcg, only: prn - use search, only: binary_search - use string, only: to_str - use tally, only: score_analog_tally, score_tracklength_tally, & - score_surface_current + use interpolation, only: interpolate_tab1 + use material_header, only: Material + use mesh, only: get_mesh_indices + use output, only: write_message + use particle_header, only: LocalCoord + use particle_restart_write, only: write_particle_restart + use random_lcg, only: prn + use search, only: binary_search + use string, only: to_str + use tally, only: score_analog_tally, & + score_tracklength_tally, & + score_surface_current implicit none @@ -54,6 +56,7 @@ contains ! Particle couldn't be located if (.not. found_cell) then + call write_particle_restart() message = "Could not locate particle " // trim(to_str(p % id)) call fatal_error() end if @@ -262,6 +265,7 @@ contains ! Check to make sure that a nuclide was sampled if (i > mat % n_nuclides) then + call write_particle_restart() message = "Did not sample any nuclide during collision." call fatal_error() end if @@ -475,6 +479,7 @@ contains ! Check to make sure inelastic scattering reaction sampled if (i > nuc % n_reaction) then + call write_particle_restart() message = "Did not sample any reaction for nuclide " // & trim(nuc % name) // " on material " // & trim(to_str(mat % id)) @@ -905,6 +910,7 @@ contains ! Determine indices on ufs mesh for current location call get_mesh_indices(ufs_mesh, p % coord0 % xyz, ijk, in_mesh) if (.not. in_mesh) then + call write_particle_restart() message = "Source site outside UFS mesh!" call fatal_error() end if @@ -997,6 +1003,7 @@ contains ! check for large number of resamples n_sample = n_sample + 1 if (n_sample == MAX_SAMPLE) then + call write_particle_restart() message = "Resampled energy distribution maximum number of " // & "times for nuclide " // nuc % name call fatal_error() @@ -1023,6 +1030,7 @@ contains ! check for large number of resamples n_sample = n_sample + 1 if (n_sample == MAX_SAMPLE) then + call write_particle_restart() message = "Resampled energy distribution maximum number of " // & "times for nuclide " // nuc % name call fatal_error() @@ -1233,6 +1241,7 @@ contains mu = mu0 + (sqrt(max(ZERO, p0*p0 + 2*frac*(xi - c_k))) - p0)/frac end if else + call write_particle_restart() message = "Unknown interpolation type: " // trim(to_str(interp)) call fatal_error() end if @@ -1244,6 +1253,7 @@ contains if (abs(mu) > ONE) mu = sign(ONE,mu) else + call write_particle_restart() message = "Unknown angular distribution type: " // trim(to_str(type)) call fatal_error() end if @@ -1394,6 +1404,7 @@ contains NE = int(edist % data(2 + 2*NR)) NET = int(edist % data(3 + 2*NR + NE)) if (NR > 0) then + call write_particle_restart() message = "Multiple interpolation regions not supported while & &attempting to sample equiprobable energy bins." call fatal_error() @@ -1459,6 +1470,7 @@ contains NR = int(edist % data(1)) NE = int(edist % data(2 + 2*NR)) if (NR > 0) then + call write_particle_restart() message = "Multiple interpolation regions not supported while & &attempting to sample continuous tabular distribution." call fatal_error() @@ -1518,6 +1530,7 @@ contains if (ND > 0) then ! discrete lines present + call write_particle_restart() message = "Discrete lines in continuous tabular distributed not & &yet supported" call fatal_error() @@ -1559,6 +1572,7 @@ contains 2*frac*(r1 - c_k))) - p_l_k)/frac end if else + call write_particle_restart() message = "Unknown interpolation type: " // trim(to_str(INTT)) call fatal_error() end if @@ -1600,6 +1614,7 @@ contains ! check for large number of rejections n_sample = n_sample + 1 if (n_sample == MAX_SAMPLE) then + call write_particle_restart() message = "Too many rejections on Maxwell fission spectrum." call fatal_error() end if @@ -1632,6 +1647,7 @@ contains ! check for large number of rejections n_sample = n_sample + 1 if (n_sample == MAX_SAMPLE) then + call write_particle_restart() message = "Too many rejections on evaporation spectrum." call fatal_error() end if @@ -1673,6 +1689,7 @@ contains ! check for large number of rejections n_sample = n_sample + 1 if (n_sample == MAX_SAMPLE) then + call write_particle_restart() message = "Too many rejections on Watt spectrum." call fatal_error() end if @@ -1683,6 +1700,7 @@ contains ! KALBACH-MANN CORRELATED SCATTERING if (.not. present(mu_out)) then + call write_particle_restart() message = "Law 44 called without giving mu_out as argument." call fatal_error() end if @@ -1691,6 +1709,7 @@ contains NR = int(edist % data(1)) NE = int(edist % data(2 + 2*NR)) if (NR > 0) then + call write_particle_restart() message = "Multiple interpolation regions not supported while & &attempting to sample Kalbach-Mann distribution." call fatal_error() @@ -1751,6 +1770,7 @@ contains if (ND > 0) then ! discrete lines present + call write_particle_restart() message = "Discrete lines in continuous tabular distributed not & &yet supported" call fatal_error() @@ -1806,6 +1826,7 @@ contains KM_R = R_k + (R_k1 - R_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k) KM_A = A_k + (A_k1 - A_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k) else + call write_particle_restart() message = "Unknown interpolation type: " // trim(to_str(INTT)) call fatal_error() end if @@ -1832,6 +1853,7 @@ contains ! CORRELATED ENERGY AND ANGLE DISTRIBUTION if (.not. present(mu_out)) then + call write_particle_restart() message = "Law 44 called without giving mu_out as argument." call fatal_error() end if @@ -1840,6 +1862,7 @@ contains NR = int(edist % data(1)) NE = int(edist % data(2 + 2*NR)) if (NR > 0) then + call write_particle_restart() message = "Multiple interpolation regions not supported while & &attempting to sample correlated energy-angle distribution." call fatal_error() @@ -1900,6 +1923,7 @@ contains if (ND > 0) then ! discrete lines present + call write_particle_restart() message = "Discrete lines in continuous tabular distributed not & &yet supported" call fatal_error() @@ -1942,6 +1966,7 @@ contains 2*frac*(r1 - c_k))) - p_l_k)/frac end if else + call write_particle_restart() message = "Unknown interpolation type: " // trim(to_str(INTT)) call fatal_error() end if @@ -2005,6 +2030,7 @@ contains mu_out = mu_k + (sqrt(p_k*p_k + 2*frac*(r3 - c_k))-p_k)/frac end if else + call write_particle_restart() message = "Unknown interpolation type: " // trim(to_str(JJ)) call fatal_error() end if From 19e9c28bb69a766c6d61be2b462cbaaa6d301580 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Mon, 8 Apr 2013 05:41:15 -0700 Subject: [PATCH 08/21] fixed #endif statement --- src/particle_restart_write.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/particle_restart_write.F90 b/src/particle_restart_write.F90 index d9388ca76..da6a4067d 100644 --- a/src/particle_restart_write.F90 +++ b/src/particle_restart_write.F90 @@ -16,7 +16,7 @@ module particle_restart_write #ifdef HDF5 integer(HID_T) :: hdf5_particle_file -#endif HDF5 +#endif contains From 72e83587b6b05af867da422b0a8d7ebaf11ce038 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Mon, 8 Apr 2013 05:46:32 -0700 Subject: [PATCH 09/21] made current_work comment more descriptive --- src/global.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/global.F90 b/src/global.F90 index 12c88fc58..33dd9b759 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -166,7 +166,7 @@ module global integer(8) :: bank_last ! index of last particle in bank integer(8) :: work ! number of particles per processor integer(8) :: maxwork ! maximum number of particles per processor - integer(8) :: current_work ! current work + integer(8) :: current_work ! index in source bank of current history simulated ! Temporary k-effective values real(8), allocatable :: k_batch(:) ! batch estimates of k From c27d3fbf481b7d5f87e8b1710ef0d8bd8ae876c6 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Mon, 8 Apr 2013 05:55:21 -0700 Subject: [PATCH 10/21] fixed typo in for particle restart option --- src/initialize.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index 33fb953d5..aa418cbe7 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -306,7 +306,7 @@ contains case ('-eps_tol', '-ksp_gmres_restart') ! Handle options that would be based to PETSC i = i + 1 - case ('-partcle','--particle') + case ('-particle','--particle') ! Read in path for particle restart i = i + 1 path_particle_restart = argv(i) From d0919de1b82c970941e16cac2f48e44ce930aa01 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Mon, 8 Apr 2013 06:08:25 -0700 Subject: [PATCH 11/21] got rid of previous loop counter i completely --- src/eigenvalue.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 9c0678f7a..660a1119f 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -38,8 +38,6 @@ contains subroutine run_eigenvalue() - integer(8) :: i ! index over individual particles - if (master) call header("K EIGENVALUE SIMULATION", level=1) ! Allocate particle From 34bb27bf7dc96093f154c54863fa66c7c59d6801 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Mon, 8 Apr 2013 08:39:03 -0700 Subject: [PATCH 12/21] added binary writer and reader for particle restart --- src/constants.F90 | 15 ++++---- src/particle_restart.F90 | 63 +++++++++++++++++++++++++++++----- src/particle_restart_write.F90 | 59 ++++++++++++++++++++++++++++++- 3 files changed, 120 insertions(+), 17 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 156f0046a..c4cbf29d0 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -350,13 +350,14 @@ module constants MODE_PARTICLE = 5 ! Particle restart mode ! Unit numbers - integer, parameter :: UNIT_SUMMARY = 11 ! unit # for writing summary file - integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file - integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file - integer, parameter :: UNIT_XS = 14 ! unit # for writing xs summary file - integer, parameter :: UNIT_SOURCE = 15 ! unit # for writing source file - integer, parameter :: UNIT_STATE = 16 ! unit # for writing state point - integer, parameter :: CMFD_BALANCE = 17 ! unit # for writing cmfd balance file + integer, parameter :: UNIT_SUMMARY = 11 ! unit # for writing summary file + integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file + integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file + integer, parameter :: UNIT_XS = 14 ! unit # for writing xs summary file + integer, parameter :: UNIT_SOURCE = 15 ! unit # for writing source file + integer, parameter :: UNIT_STATE = 16 ! unit # for writing state point + integer, parameter :: CMFD_BALANCE = 17 ! unit # for writing cmfd balance file + integer, parameter :: UNIT_PARTICLE = 18 ! unit # for writing particle restart !============================================================================= ! CMFD CONSTANTS diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index ac2730b8f..93e1d8084 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -30,12 +30,14 @@ module particle_restart contains +#ifdef HDF5 + !=============================================================================== -! READ_PARTICLE_RESTART +! READ_HDF5_PARTICLE_RESTART !=============================================================================== subroutine read_hdf5_particle_restart() -#ifdef HDF5 + integer(HSIZE_T) :: dims1(1) ! write meessage @@ -68,15 +70,52 @@ contains ! close hdf5 file call h5fclose_f(hdf5_particle_file, hdf5_err) -#endif + end subroutine read_hdf5_particle_restart +#endif + +!=============================================================================== +! READ_BINARY_PARTICLE_RESTART +!=============================================================================== + + subroutine read_binary_particle_restart() + + ! write meessage + message = "Loading particle restart file " // trim(path_particle_restart) & + // "..." + call write_message(1) + + ! open file + open(UNIT=UNIT_PARTICLE, FILE=path_particle_restart, STATUS='old', & + ACCESS='stream') + + ! read data from file + read(UNIT_PARTICLE) current_batch + read(UNIT_PARTICLE) gen_per_batch + read(UNIT_PARTICLE) current_gen + read(UNIT_PARTICLE) n_particles + read(UNIT_PARTICLE) p % id + read(UNIT_PARTICLE) p % wgt + read(UNIT_PARTICLE) p % E + read(UNIT_PARTICLE) p % coord % xyz + read(UNIT_PARTICLE) p % coord % uvw + + ! set particle last attributes + p % last_wgt = p % wgt + p % last_xyz = p % coord % xyz + p % last_E = p % E + + ! close hdf5 file + close(UNIT_PARTICLE) + + end subroutine read_binary_particle_restart + !=============================================================================== ! RUN_PARTICLE_RESTART !=============================================================================== subroutine run_particle_restart() -#ifdef HDF5 integer(8) :: particle_seed @@ -85,7 +124,11 @@ contains call initialize_particle() ! read in the restart information +#ifdef HDF5 call read_hdf5_particle_restart() +#else + call read_binary_particle_restart() +#endif ! set all tallies to 0 for now (just tracking errors) n_tallies = 0 @@ -97,12 +140,14 @@ contains ! transport neutron call transport() - print *, 'WEIGHT:', p % wgt - print *, 'ENERGY:', p % E - print *, 'LOCATION:', p % coord % xyz - print *, 'ANGLE:', p % coord % uvw -#endif HDF5 + ! write output if particle made it + write(ou,*) 'Particle Successfully Transport:' + write(ou,*) 'WEIGHT:', p % wgt + write(ou,*) 'ENERGY:', p % E + write(ou,*) 'LOCATION:', p % coord % xyz + write(ou,*) 'ANGLE:', p % coord % uvw + end subroutine run_particle_restart end module particle_restart diff --git a/src/particle_restart_write.F90 b/src/particle_restart_write.F90 index da6a4067d..4720c9dd4 100644 --- a/src/particle_restart_write.F90 +++ b/src/particle_restart_write.F90 @@ -25,7 +25,27 @@ contains !=============================================================================== subroutine write_particle_restart() + + ! Dont write another restart file if in particle restart mode + if (run_mode == MODE_PARTICLE) return + + ! write out binary or HDF5 file #ifdef HDF5 + call write_particle_restart_hdf5() +#else + call write_particle_restart_binary() +#endif + + end subroutine write_particle_restart + +#ifdef HDF5 + +!=============================================================================== +! WRITE_PARTICLE_RESTART_HDF5 +!=============================================================================== + + subroutine write_particle_restart_hdf5() + character(MAX_FILE_LEN) :: filename integer(HSIZE_T) :: dims1(1) type(Bank), pointer :: src => null() @@ -55,7 +75,44 @@ contains ! close hdf5 file call h5fclose_f(hdf5_particle_file, hdf5_err) + + end subroutine write_particle_restart_hdf5 + #endif - end subroutine write_particle_restart + +!=============================================================================== +! WRITE_PARTICLE_RESTART_BINARY +!=============================================================================== + + subroutine write_particle_restart_binary() + + character(MAX_FILE_LEN) :: filename + type(Bank), pointer :: src => null() + + ! set up file name + filename = 'particle_'//trim(to_str(rank))//'.binary' + + ! create hdf5 file + open(UNIT=UNIT_PARTICLE, FILE=filename, STATUS='replace', & + ACCESS='stream') + + ! get information about source particle + src => source_bank(current_work) + + ! write data to file + write(UNIT_PARTICLE) current_batch + write(UNIT_PARTICLE) gen_per_batch + write(UNIT_PARTICLE) current_gen + write(UNIT_PARTICLE) n_particles + write(UNIT_PARTICLE) p % id + write(UNIT_PARTICLE) src % wgt + write(UNIT_PARTICLE) src % E + write(UNIT_PARTICLE) src % xyz + write(UNIT_PARTICLE) src % uvw + + ! close hdf5 file + close(UNIT_PARTICLE) + + end subroutine write_particle_restart_binary end module particle_restart_write From ac1397b3d14f825703b88d20f99bd4cb55d678fc Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Mon, 8 Apr 2013 08:40:54 -0700 Subject: [PATCH 13/21] only 1 particle can be restarted at a time and should run from master --- src/main.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main.F90 b/src/main.F90 index 41a7481da..e0182902a 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -26,7 +26,7 @@ program main ! For tallies-only mode, we just skip straight to finalize_run to write out ! the tally results case (MODE_PARTICLE) - call run_particle_restart() + if(master) call run_particle_restart() end select ! finalize run From ceb5730030e1cc8837465227928825e9972e8baa Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Mon, 8 Apr 2013 08:44:39 -0700 Subject: [PATCH 14/21] added test for particle restart --- tests/test_particle_restart/geometry.xml | 12 +++++++ tests/test_particle_restart/materials.xml | 9 ++++++ tests/test_particle_restart/settings.xml | 16 ++++++++++ .../test_particle_restart.py | 31 +++++++++++++++++++ 4 files changed, 68 insertions(+) create mode 100644 tests/test_particle_restart/geometry.xml create mode 100644 tests/test_particle_restart/materials.xml create mode 100644 tests/test_particle_restart/settings.xml create mode 100644 tests/test_particle_restart/test_particle_restart.py diff --git a/tests/test_particle_restart/geometry.xml b/tests/test_particle_restart/geometry.xml new file mode 100644 index 000000000..3dc9f24c1 --- /dev/null +++ b/tests/test_particle_restart/geometry.xml @@ -0,0 +1,12 @@ + + + + + + + + + + + + diff --git a/tests/test_particle_restart/materials.xml b/tests/test_particle_restart/materials.xml new file mode 100644 index 000000000..1f85510e0 --- /dev/null +++ b/tests/test_particle_restart/materials.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/test_particle_restart/settings.xml b/tests/test_particle_restart/settings.xml new file mode 100644 index 000000000..f278bb5b0 --- /dev/null +++ b/tests/test_particle_restart/settings.xml @@ -0,0 +1,16 @@ + + + + + 10 + 5 + 1000 + + + + + -20 -20 -5 20 20 5 + + + + diff --git a/tests/test_particle_restart/test_particle_restart.py b/tests/test_particle_restart/test_particle_restart.py new file mode 100644 index 000000000..5c82f1824 --- /dev/null +++ b/tests/test_particle_restart/test_particle_restart.py @@ -0,0 +1,31 @@ +#!/usr/bin/env python + +import os +from subprocess import Popen, STDOUT, PIPE + +pwd = os.path.dirname(__file__) + +def setup(): + os.putenv('PWD', pwd) + os.chdir(pwd) + +def test_run(): + proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE) + returncode = proc.wait() + print(proc.communicate()[0]) + assert returncode == 0 + +def test_run_restart(): + proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE) + returncode = proc.wait() + print(proc.communicate()[0]) + assert returncode == 0 + +def test_created_restart(): + assert os.path.exists(pwd + '/particle_0.binary') + +def teardown(): + output = [pwd + '/particle_0.binary'] + for f in output: + if os.path.exists(f): + os.remove(f) From d3afc2f64cc7d43b26dcab6fb1a9eafebacc3ec9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 8 Apr 2013 20:01:46 -0400 Subject: [PATCH 15/21] Fixed another #endif statement. --- src/particle_restart.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 93e1d8084..07c7f5475 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -22,7 +22,7 @@ module particle_restart #ifdef HDF5 integer(HID_T) :: hdf5_particle_file -#endif HDF5 +#endif ! Short names for output and error units integer :: ou = OUTPUT_UNIT From 878697dfc9f2ff49454f36c6517e3f1830f2a344 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 8 Apr 2013 20:04:23 -0400 Subject: [PATCH 16/21] Display correct surface id on lost particle error. --- src/geometry.F90 | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 3ab85ba72..685971288 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -247,22 +247,24 @@ contains integer, intent(in) :: last_cell ! last cell particle was in - real(8) :: x ! x-x0 for sphere - real(8) :: y ! y-y0 for sphere - real(8) :: z ! z-z0 for sphere - real(8) :: R ! radius of sphere - real(8) :: u ! x-component of direction - real(8) :: v ! y-component of direction - real(8) :: w ! z-component of direction - real(8) :: n1 ! x-component of surface normal - real(8) :: n2 ! y-component of surface normal - real(8) :: n3 ! z-component of surface normal - real(8) :: dot_prod ! dot product of direction and normal - real(8) :: norm ! "norm" of surface normal - logical :: found ! particle found in universe? + real(8) :: x ! x-x0 for sphere + real(8) :: y ! y-y0 for sphere + real(8) :: z ! z-z0 for sphere + real(8) :: R ! radius of sphere + real(8) :: u ! x-component of direction + real(8) :: v ! y-component of direction + real(8) :: w ! z-component of direction + real(8) :: n1 ! x-component of surface normal + real(8) :: n2 ! y-component of surface normal + real(8) :: n3 ! z-component of surface normal + real(8) :: dot_prod ! dot product of direction and normal + real(8) :: norm ! "norm" of surface normal + integer :: i_surface ! index in surfaces + logical :: found ! particle found in universe? type(Surface), pointer :: surf => null() - surf => surfaces(abs(p % surface)) + i_surface = abs(p % surface) + surf => surfaces(i_surface) if (verbosity >= 10 .or. trace) then message = " Crossing surface " // trim(to_str(surf % id)) call write_message() @@ -521,7 +523,7 @@ contains if (.not. found) then call write_particle_restart() message = "After particle " // trim(to_str(p % id)) // " crossed surface " & - // trim(to_str(surfaces(abs(p%surface)) % id)) // " it could not be & + // trim(to_str(surfaces(i_surface) % id)) // " it could not be & &located in any cell and it did not leak." call fatal_error() end if From 0b8312a0c9f481c4485061f851f8e0c514f1ae56 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 8 Apr 2013 20:17:10 -0400 Subject: [PATCH 17/21] Add space between if and ( --- src/main.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main.F90 b/src/main.F90 index e0182902a..5e5ea89ba 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -26,7 +26,7 @@ program main ! For tallies-only mode, we just skip straight to finalize_run to write out ! the tally results case (MODE_PARTICLE) - if(master) call run_particle_restart() + if (master) call run_particle_restart() end select ! finalize run From 86594ee03ae93a80c42beb2add79a72c50050ae9 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 9 Apr 2013 05:20:54 -0700 Subject: [PATCH 18/21] -s command option can now be used for particle restart --- src/initialize.F90 | 2 +- src/output.F90 | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index bdee6b780..a1518cd21 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -306,7 +306,7 @@ contains case ('-eps_tol', '-ksp_gmres_restart') ! Handle options that would be based to PETSC i = i + 1 - case ('-particle','--particle') + case ('-s','-particle','--particle') ! Read in path for particle restart i = i + 1 path_particle_restart = argv(i) diff --git a/src/output.F90 b/src/output.F90 index ec646ccec..b5b1d1168 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -167,6 +167,7 @@ contains write(OUTPUT_UNIT,*) 'Options:' write(OUTPUT_UNIT,*) ' -p, --plot Run in plotting mode' write(OUTPUT_UNIT,*) ' -r, --restart Restart a previous run' + write(OUTPUT_UNIT,*) ' -s, --particle Run a single particle history' write(OUTPUT_UNIT,*) ' -t, --tallies Write tally results from state point' write(OUTPUT_UNIT,*) ' -v, --version Show version information' write(OUTPUT_UNIT,*) ' -?, --help Show this message' From 8a83d1a4d426a82505fc93ad032e2a87fbba2117 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 9 Apr 2013 05:26:14 -0700 Subject: [PATCH 19/21] particle restart test uses -s option, return codes not equal to 0 --- tests/test_particle_restart/test_particle_restart.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/test_particle_restart/test_particle_restart.py b/tests/test_particle_restart/test_particle_restart.py index 5c82f1824..a69065f1b 100644 --- a/tests/test_particle_restart/test_particle_restart.py +++ b/tests/test_particle_restart/test_particle_restart.py @@ -13,13 +13,14 @@ def test_run(): proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE) returncode = proc.wait() print(proc.communicate()[0]) - assert returncode == 0 + assert returncode != 0 def test_run_restart(): - proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE) + proc = Popen([pwd + '/../../src/openmc -s particle_0.binary'], + stderr=STDOUT, stdout=PIPE) returncode = proc.wait() print(proc.communicate()[0]) - assert returncode == 0 + assert returncode != 0 def test_created_restart(): assert os.path.exists(pwd + '/particle_0.binary') From 21cbbe823f7c09c582f8024f12dfc098c2847124 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 9 Apr 2013 16:48:34 -0400 Subject: [PATCH 20/21] Added -s, --particle option to man page. --- man/man1/openmc.1 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/man/man1/openmc.1 b/man/man1/openmc.1 index fb7adb518..b8a9d9c44 100644 --- a/man/man1/openmc.1 +++ b/man/man1/openmc.1 @@ -19,6 +19,10 @@ Run in plotting mode .BI \-r " statepoint" "\fR,\fP \-\-restart" " statepoint" Restart a previous run using data stored in \fIstatepoint\fP. .TP +.BI \-s " particleFile" "\fR,\fP \-\-particle" " particleFile" +Run a single particle history from a particle restart file named +\fIparticleFile\fP. +.TP .B "\-v\fR, \fP\-\-version" Show version information. .TP From f83a2dd996ad8f1f0807e4ef877df52d62a4ce45 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 9 Apr 2013 17:19:10 -0400 Subject: [PATCH 21/21] Change particle restart test to check for output on stderr rather than checking the return status. --- .../test_particle_restart.py | 20 +++++++++---------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/tests/test_particle_restart/test_particle_restart.py b/tests/test_particle_restart/test_particle_restart.py index a69065f1b..b1f9a0108 100644 --- a/tests/test_particle_restart/test_particle_restart.py +++ b/tests/test_particle_restart/test_particle_restart.py @@ -10,21 +10,19 @@ def setup(): os.chdir(pwd) def test_run(): - proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE) - returncode = proc.wait() - print(proc.communicate()[0]) - assert returncode != 0 - -def test_run_restart(): - proc = Popen([pwd + '/../../src/openmc -s particle_0.binary'], - stderr=STDOUT, stdout=PIPE) - returncode = proc.wait() - print(proc.communicate()[0]) - assert returncode != 0 + proc = Popen([pwd + '/../../src/openmc'], stderr=PIPE, stdout=PIPE) + stdout, stderr = proc.communicate() + assert stderr != '' def test_created_restart(): assert os.path.exists(pwd + '/particle_0.binary') +def test_run_restart(): + proc = Popen([pwd + '/../../src/openmc -s particle_0.binary'], + stderr=PIPE, stdout=PIPE, shell=True) + stdout, stderr = proc.communicate() + assert stderr != '' + def teardown(): output = [pwd + '/particle_0.binary'] for f in output: