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
diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES
index 61c31dd90..5e2a96a09 100644
--- a/src/DEPENDENCIES
+++ b/src/DEPENDENCIES
@@ -154,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
@@ -239,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
@@ -264,6 +266,22 @@ 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: 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
physics.o: cross_section.o
@@ -278,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 d53f0a543..c9612e6af 100644
--- a/src/OBJECTS
+++ b/src/OBJECTS
@@ -40,6 +40,8 @@ mesh_header.o \
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/constants.F90 b/src/constants.F90
index 1a3d330ff..c4cbf29d0 100644
--- a/src/constants.F90
+++ b/src/constants.F90
@@ -346,16 +346,18 @@ 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
- 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/eigenvalue.F90 b/src/eigenvalue.F90
index 8357c89f3..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
@@ -74,10 +72,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/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/geometry.F90 b/src/geometry.F90
index 64912cb95..685971288 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
@@ -246,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()
@@ -303,6 +306,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 +440,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 +461,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,8 +521,9 @@ 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 &
+ // 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
@@ -608,6 +615,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 +638,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/global.F90 b/src/global.F90
index 005426211..5042ac87e 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(8) :: current_work ! index in source bank of current history simulated
! Temporary k-effective values
real(8), allocatable :: k_batch(:) ! batch estimates of k
@@ -253,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
@@ -274,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 83946165d..a1518cd21 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 ('-s','-particle','--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..5e5ea89ba 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)
+ if (master) call run_particle_restart()
end select
! finalize run
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'
diff --git a/src/particle_restart.F90 b/src/particle_restart.F90
new file mode 100644
index 000000000..07c7f5475
--- /dev/null
+++ b/src/particle_restart.F90
@@ -0,0 +1,153 @@
+module particle_restart
+
+ use, intrinsic :: ISO_FORTRAN_ENV
+
+ use bank_header, only: Bank
+ use constants
+ 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_interface
+#endif
+
+ implicit none
+ private
+ public :: run_particle_restart
+
+#ifdef HDF5
+ integer(HID_T) :: hdf5_particle_file
+#endif
+
+ ! Short names for output and error units
+ integer :: ou = OUTPUT_UNIT
+ integer :: eu = ERROR_UNIT
+
+contains
+
+#ifdef HDF5
+
+!===============================================================================
+! READ_HDF5_PARTICLE_RESTART
+!===============================================================================
+
+ subroutine read_hdf5_particle_restart()
+
+ 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)
+
+ ! 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_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()
+
+ integer(8) :: particle_seed
+
+ ! initialize the particle to be tracked
+ allocate(p)
+ 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
+
+ ! 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()
+
+ ! 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
new file mode 100644
index 000000000..4720c9dd4
--- /dev/null
+++ b/src/particle_restart_write.F90
@@ -0,0 +1,118 @@
+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
+
+contains
+
+!===============================================================================
+! WRITE_PARTICLE_RESTART
+!===============================================================================
+
+ 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()
+
+ ! 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)
+
+ end subroutine write_particle_restart_hdf5
+
+#endif
+
+!===============================================================================
+! 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
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
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..b1f9a0108
--- /dev/null
+++ b/tests/test_particle_restart/test_particle_restart.py
@@ -0,0 +1,30 @@
+#!/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=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:
+ if os.path.exists(f):
+ os.remove(f)