Move particle restart writing / lost particle handling to particle_header

This commit is contained in:
Paul Romano 2017-12-11 14:02:45 +07:00
parent a149ef42a0
commit d4f366cbd4
8 changed files with 132 additions and 147 deletions

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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