From d4f366cbd40d05f0981ef890ed0e61ae2b3d91f4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 11 Dec 2017 14:02:45 +0700 Subject: [PATCH] Move particle restart writing / lost particle handling to particle_header --- CMakeLists.txt | 1 - src/api.F90 | 1 + src/geometry.F90 | 62 ++++------------ src/particle_header.F90 | 125 +++++++++++++++++++++++++++++---- src/particle_restart_write.F90 | 74 ------------------- src/physics.F90 | 11 ++- src/physics_mg.F90 | 3 +- src/simulation_header.F90 | 2 +- 8 files changed, 132 insertions(+), 147 deletions(-) delete mode 100644 src/particle_restart_write.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 8ca62b3e5..c5150745b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -377,7 +377,6 @@ set(LIBOPENMC_FORTRAN_SRC src/output.F90 src/particle_header.F90 src/particle_restart.F90 - src/particle_restart_write.F90 src/physics_common.F90 src/physics.F90 src/physics_mg.F90 diff --git a/src/api.F90 b/src/api.F90 index 4538bbf3e..326b755b6 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -112,6 +112,7 @@ contains legendre_to_tabular = .true. legendre_to_tabular_points = 33 n_batch_interval = 1 + n_lost_particles = 0 n_particles = 0 n_source_points = 0 n_state_points = 0 diff --git a/src/geometry.F90 b/src/geometry.F90 index 1cc27e09e..1d99273de 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -5,7 +5,6 @@ module geometry use geometry_header use output, only: write_message use particle_header, only: LocalCoord, Particle - use particle_restart_write, only: write_particle_restart use simulation_header use settings use surface_header @@ -357,7 +356,7 @@ contains else ! Particle is outside the lattice. if (lat % outer == NO_OUTER_UNIVERSE) then - call handle_lost_particle(p, "Particle " // trim(to_str(p %id)) & + call p % mark_as_lost("Particle " // trim(to_str(p %id)) & // " is outside lattice " // trim(to_str(lat % id)) & // " but the lattice has no defined outer universe.") return @@ -440,8 +439,8 @@ contains ! Do not handle reflective boundary conditions on lower universes if (p % n_coord /= 1) then - call handle_lost_particle(p, "Cannot reflect particle " & - &// trim(to_str(p % id)) // " off surface in a lower universe.") + call p % mark_as_lost("Cannot reflect particle " & + // trim(to_str(p % id)) // " off surface in a lower universe.") return end if @@ -477,8 +476,8 @@ contains p % n_coord = 1 call find_cell(p, found) if (.not. found) then - call handle_lost_particle(p, "Couldn't find particle after reflecting& - & from surface " // trim(to_str(surf%id)) // ".") + call p % mark_as_lost("Couldn't find particle after reflecting& + & from surface " // trim(to_str(surf % id)) // ".") return end if @@ -497,7 +496,7 @@ contains ! Do not handle periodic boundary conditions on lower universes if (p % n_coord /= 1) then - call handle_lost_particle(p, "Cannot transfer particle " & + call p % mark_as_lost("Cannot transfer particle " & // trim(to_str(p % id)) // " across surface in a lower universe.& & Boundary conditions must be applied to universe 0.") return @@ -584,8 +583,8 @@ contains p % n_coord = 1 call find_cell(p, found) if (.not. found) then - call handle_lost_particle(p, "Couldn't find particle after hitting & - &periodic boundary on surface " // trim(to_str(surf%id)) // ".") + call p % mark_as_lost("Couldn't find particle after hitting & + &periodic boundary on surface " // trim(to_str(surf % id)) // ".") return end if @@ -641,8 +640,8 @@ contains ! undefined region in the geometry. if (.not. found) then - call handle_lost_particle(p, "After particle " // trim(to_str(p % id)) & - // " crossed surface " // trim(to_str(surf%id)) & + call p % mark_as_lost("After particle " // trim(to_str(p % id)) & + // " crossed surface " // trim(to_str(surf % id)) & // " it could not be located in any cell and it did not leak.") return end if @@ -690,7 +689,7 @@ contains call find_cell(p, found) if (.not. found) then if (p % alive) then ! Particle may have been killed in find_cell - call handle_lost_particle(p, "Could not locate particle " & + call p % mark_as_lost("Could not locate particle " & // trim(to_str(p % id)) // " after crossing a lattice boundary.") return end if @@ -713,9 +712,8 @@ contains ! Search for particle call find_cell(p, found) if (.not. found) then - call handle_lost_particle(p, "Could not locate particle " & - // trim(to_str(p % id)) & - // " after crossing a lattice boundary.") + call p % mark_as_lost("Could not locate particle " // & + trim(to_str(p % id)) // " after crossing a lattice boundary.") return end if end if @@ -992,7 +990,7 @@ contains end select LAT_TYPE if (d_lat < ZERO) then - call handle_lost_particle(p, "Particle " // trim(to_str(p % id)) & + call p % mark_as_lost("Particle " // trim(to_str(p % id)) & //" had a negative distance to a lattice boundary. d = " & //trim(to_str(d_lat))) end if @@ -1093,38 +1091,6 @@ contains end subroutine neighbor_lists -!=============================================================================== -! HANDLE_LOST_PARTICLE -!=============================================================================== - - subroutine handle_lost_particle(p, message) - - type(Particle), intent(inout) :: p - character(*) :: message - - integer(8) :: tot_n_particles - - ! Print warning and write lost particle file - call warning(message) - call write_particle_restart(p) - - ! Increment number of lost particles - p % alive = .false. -!$omp atomic - n_lost_particles = n_lost_particles + 1 - - ! Count the total number of simulated particles (on this processor) - tot_n_particles = n_batches * gen_per_batch * work - - ! Abort the simulation if the maximum number of lost particles has been - ! reached - if (n_lost_particles >= MAX_LOST_PARTICLES .and. & - n_lost_particles >= REL_MAX_LOST_PARTICLES * tot_n_particles) then - call fatal_error("Maximum number of lost particles has been reached.") - end if - - end subroutine handle_lost_particle - !=============================================================================== ! CALC_OFFSETS calculates and stores the offsets in all fill cells. This ! routine is called once upon initialization. diff --git a/src/particle_header.F90 b/src/particle_header.F90 index bda7f03b4..66b50c134 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -1,20 +1,27 @@ module particle_header - use bank_header, only: Bank - use constants, only: NEUTRON, ONE, NONE, ZERO, MAX_SECONDARY, & - MAX_DELAYED_GROUPS, ERROR_REAL - use error, only: fatal_error + use hdf5, only: HID_T + + use bank_header, only: Bank, source_bank + use constants + use error, only: fatal_error, warning use geometry_header, only: root_universe + use hdf5_interface + use settings + use simulation_header + use string, only: to_str implicit none + private + !=============================================================================== ! LOCALCOORD describes the location of a particle local to a single ! universe. When the geometry consists of nested universes, a particle will have ! a list of coordinates in each level !=============================================================================== - type LocalCoord + type, public :: LocalCoord ! Indices in various arrays for this level integer :: cell = NONE @@ -39,7 +46,7 @@ module particle_header ! geometry !=============================================================================== - type Particle + type, public :: Particle ! Basic data integer(8) :: id ! Unique ID integer :: type ! Particle type (n, p, e, etc) @@ -107,10 +114,12 @@ module particle_header type(Bank) :: secondary_bank(MAX_SECONDARY) contains - procedure :: initialize => initialize_particle - procedure :: clear => clear_particle - procedure :: initialize_from_source + procedure :: clear procedure :: create_secondary + procedure :: initialize + procedure :: initialize_from_source + procedure :: mark_as_lost + procedure :: write_restart end type Particle contains @@ -120,7 +129,7 @@ contains ! bank !=============================================================================== - subroutine initialize_particle(this) + subroutine initialize(this) class(Particle) :: this @@ -154,23 +163,22 @@ contains this % n_coord = 1 this % last_n_coord = 1 - end subroutine initialize_particle + end subroutine initialize !=============================================================================== ! CLEAR_PARTICLE resets all coordinate levels for the particle !=============================================================================== - subroutine clear_particle(this) - + subroutine clear(this) class(Particle) :: this + integer :: i ! remove any coordinate levels do i = 1, MAX_COORD call this % coord(i) % reset() end do - - end subroutine clear_particle + end subroutine clear !=============================================================================== ! RESET_COORD clears data from a single coordinate level @@ -255,4 +263,91 @@ contains end subroutine create_secondary +!=============================================================================== +! MARK_AS_LOST +!=============================================================================== + + subroutine mark_as_lost(this, message) + class(Particle), intent(inout) :: this + character(*) :: message + + integer(8) :: tot_n_particles + + ! Print warning and write lost particle file + call warning(message) + call this % write_restart() + + ! Increment number of lost particles + this % alive = .false. +!$omp atomic + n_lost_particles = n_lost_particles + 1 + + ! Count the total number of simulated particles (on this processor) + tot_n_particles = current_batch * gen_per_batch * work + + ! Abort the simulation if the maximum number of lost particles has been + ! reached + if (n_lost_particles >= MAX_LOST_PARTICLES .and. & + n_lost_particles >= REL_MAX_LOST_PARTICLES * tot_n_particles) then + call fatal_error("Maximum number of lost particles has been reached.") + end if + + end subroutine mark_as_lost + +!=============================================================================== +! WRITE_RESTART creates a particle restart file +!=============================================================================== + + subroutine write_restart(this) + class(Particle), intent(in) :: this + + integer(HID_T) :: file_id + character(MAX_FILE_LEN) :: filename + + ! 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(this % id)) // '.h5' + +!$omp critical (WriteParticleRestart) + ! Create file + file_id = file_create(filename) + + associate (src => source_bank(current_work)) + ! Write filetype and version info + call write_attribute(file_id, 'filetype', 'particle restart') + call write_attribute(file_id, 'version', VERSION_PARTICLE_RESTART) + call write_attribute(file_id, "openmc_version", VERSION) +#ifdef GIT_SHA1 + call write_attribute(file_id, "git_sha1", GIT_SHA1) +#endif + + ! Write data to file + call write_dataset(file_id, 'current_batch', current_batch) + call write_dataset(file_id, 'generations_per_batch', gen_per_batch) + call write_dataset(file_id, 'current_generation', current_gen) + call write_dataset(file_id, 'n_particles', n_particles) + select case(run_mode) + case (MODE_FIXEDSOURCE) + call write_dataset(file_id, 'run_mode', 'fixed source') + case (MODE_EIGENVALUE) + call write_dataset(file_id, 'run_mode', 'eigenvalue') + case (MODE_PARTICLE) + call write_dataset(file_id, 'run_mode', 'particle restart') + end select + call write_dataset(file_id, 'id', this % 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) + end associate + + ! Close file + call file_close(file_id) +!$omp end critical (WriteParticleRestart) + + end subroutine write_restart + end module particle_header diff --git a/src/particle_restart_write.F90 b/src/particle_restart_write.F90 deleted file mode 100644 index 5addf6056..000000000 --- a/src/particle_restart_write.F90 +++ /dev/null @@ -1,74 +0,0 @@ -module particle_restart_write - - use bank_header, only: Bank, source_bank - use hdf5_interface - use particle_header, only: Particle - use settings - use simulation_header - use string, only: to_str - - use hdf5 - - implicit none - private - public :: write_particle_restart - -contains - -!=============================================================================== -! WRITE_PARTICLE_RESTART is the main routine that writes out the particle file -!=============================================================================== - - subroutine write_particle_restart(p) - type(Particle), intent(in) :: p - - integer(HID_T) :: file_id - character(MAX_FILE_LEN) :: filename - - ! 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)) // '.h5' - -!$omp critical (WriteParticleRestart) - ! Create file - file_id = file_create(filename) - - associate (src => source_bank(current_work)) - ! Write filetype and version info - call write_attribute(file_id, 'filetype', 'particle restart') - call write_attribute(file_id, 'version', VERSION_PARTICLE_RESTART) - call write_attribute(file_id, "openmc_version", VERSION) -#ifdef GIT_SHA1 - call write_attribute(file_id, "git_sha1", GIT_SHA1) -#endif - - ! Write data to file - call write_dataset(file_id, 'current_batch', current_batch) - call write_dataset(file_id, 'generations_per_batch', gen_per_batch) - call write_dataset(file_id, 'current_generation', current_gen) - call write_dataset(file_id, 'n_particles', n_particles) - select case(run_mode) - case (MODE_FIXEDSOURCE) - call write_dataset(file_id, 'run_mode', 'fixed source') - case (MODE_EIGENVALUE) - call write_dataset(file_id, 'run_mode', 'eigenvalue') - case (MODE_PARTICLE) - call write_dataset(file_id, 'run_mode', 'particle restart') - end select - 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) - end associate - - ! Close file - call file_close(file_id) -!$omp end critical (WriteParticleRestart) - - end subroutine write_particle_restart - -end module particle_restart_write diff --git a/src/physics.F90 b/src/physics.F90 index fd6a0087c..96dabd542 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -12,7 +12,6 @@ module physics use nuclide_header use output, only: write_message use particle_header, only: Particle - use particle_restart_write, only: write_particle_restart use physics_common use random_lcg, only: prn, advance_prn_seed, prn_set_stream use reaction_header, only: Reaction @@ -172,7 +171,7 @@ contains ! Check to make sure that a nuclide was sampled if (i_nuc_mat > mat % n_nuclides) then - call write_particle_restart(p) + call p % write_restart() call fatal_error("Did not sample any nuclide during collision.") end if @@ -384,7 +383,7 @@ contains ! Check to make sure inelastic scattering reaction sampled if (i > size(nuc % reactions)) then - call write_particle_restart(p) + call p % write_restart() call fatal_error("Did not sample any reaction for nuclide " & &// trim(nuc % name)) end if @@ -1098,7 +1097,7 @@ contains ! Determine indices on ufs mesh for current location call m % get_bin(p % coord(1) % xyz, mesh_bin) if (mesh_bin == NO_BIN_FOUND) then - call write_particle_restart(p) + call p % write_restart() call fatal_error("Source site outside UFS mesh!") end if @@ -1251,7 +1250,7 @@ contains ! check for large number of resamples n_sample = n_sample + 1 if (n_sample == MAX_SAMPLE) then - ! call write_particle_restart(p) + ! call p % write_restart() call fatal_error("Resampled energy distribution maximum number of " & // "times for nuclide " // nuc % name) end if @@ -1275,7 +1274,7 @@ contains ! check for large number of resamples n_sample = n_sample + 1 if (n_sample == MAX_SAMPLE) then - ! call write_particle_restart(p) + ! call p % write_restart() call fatal_error("Resampled energy distribution maximum number of " & // "times for nuclide " // nuc % name) end if diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 42094f6d0..53a28cd4d 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -13,7 +13,6 @@ module physics_mg use nuclide_header, only: material_xs use output, only: write_message use particle_header, only: Particle - use particle_restart_write, only: write_particle_restart use physics_common use random_lcg, only: prn use scattdata_header @@ -196,7 +195,7 @@ contains call m % get_bin(p % coord(1) % xyz, mesh_bin) if (mesh_bin == NO_BIN_FOUND) then - call write_particle_restart(p) + call p % write_restart() call fatal_error("Source site outside UFS mesh!") end if diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 5041f0793..825b88129 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -13,7 +13,7 @@ module simulation_header ! GEOMETRY-RELATED VARIABLES ! Number of lost particles - integer :: n_lost_particles + integer :: n_lost_particles = 0 real(8) :: log_spacing ! spacing on logarithmic grid