crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #160 from mit-crpg/particle

Browse source 
Particle Restart Files
This commit is contained in:
Paul Romano 2013-04-09 14:23:16 -07:00
commit 9a02c1558a
18 changed files with 478 additions and 62 deletions
Download patch file Download diff file Expand all files Collapse all files

4
man/man1/openmc.1
View file

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

19
src/DEPENDENCIES
View file

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

2
src/OBJECTS
View file

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

18
src/constants.F90
View file

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

6
src/eigenvalue.F90
View file

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

5
src/finalize.F90
View file

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

51
src/geometry.F90
View file

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

13
src/global.F90
View file

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

8
src/initialize.F90
View file

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

13
src/main.F90
View file

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

1
src/output.F90
View file

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

153
src/particle_restart.F90 Normal file
View file

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

118
src/particle_restart_write.F90 Normal file
View file

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

62
src/physics.F90
View file

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

12
tests/test_particle_restart/geometry.xml Normal file
View file

@ -0,0 +1,12 @@
<?xml version="1.0"?>
<geometry>
<!-- Nested Cylinders with Tracking Error -->
<surface id="1" type="z-cylinder" coeffs="0.0 0.0 10.0" />
<surface id="2" type="z-cylinder" coeffs="0.0 0.0 20.0" boundary="reflective" />
<surface id="3" type="z-plane" coeffs="5.0" boundary="reflective" />
<surface id="4" type="z-plane" coeffs="-5.0" boundary="reflective" />
<cell id="1" material="1" surfaces="-1 3 -4" />
<cell id="2" material="1" surfaces="1 -2"/>
</geometry>

9
tests/test_particle_restart/materials.xml Normal file
View file

@ -0,0 +1,9 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" xs="70c" ao="1.0" />
</material>
</materials>

16
tests/test_particle_restart/settings.xml Normal file
View file

@ -0,0 +1,16 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-20 -20 -5 20 20 5</parameters>
</space>
</source>
</settings>

30
tests/test_particle_restart/test_particle_restart.py Normal file
View file

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