From 44999762e7deedb1cb8bb7274ca33d3579a639ec Mon Sep 17 00:00:00 2001 From: Sterling Date: Wed, 12 Jun 2013 16:41:09 -0400 Subject: [PATCH] Made minor edits to implimentation. Moved particle track subroutines to the output module. Particle tracks can now be written in any run mode that calls transport. track.py should correctly handle any number of input particle tracks. --- src/DEPENDENCIES | 1 - src/output.F90 | 40 ++++++++++++++++++++++++++++++++++++++++ src/particle_restart.F90 | 26 -------------------------- src/physics.F90 | 22 ++++++++++++++++++---- src/utils/track.py | 36 +++++++++++++++++++++--------------- 5 files changed, 79 insertions(+), 46 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index b71354ccaa..fcabd989e0 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -297,7 +297,6 @@ physics.o: material_header.o physics.o: mesh.o physics.o: output.o physics.o: particle_header.o -physics.o: particle_restart.o physics.o: particle_restart_write.o physics.o: random_lcg.o physics.o: search.o diff --git a/src/output.F90 b/src/output.F90 index b8302226d7..e2114c8b48 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1945,4 +1945,44 @@ contains end function get_label +!=============================================================================== +! INITIALIZE_PARTICLE_TRACK opens a particle track output file. +! +! TODO: This subroutine needs to be modified to work with HDF5 files. It +! should also probably write a header that identifies the file as a particle +! track and maybe adds particle identifying information (batch #, etc.). +!=============================================================================== + + subroutine initialize_particle_track() + character(MAX_FILE_LEN) :: filename + + filename = trim(path_output) // 'track_' // trim(to_str(current_batch)) & + // '_' // trim(to_str(current_work)) // '.binary' + open(UNIT=UNIT_TRACK, FILE=filename, STATUS='replace', & + ACCESS='stream') + + end subroutine initialize_particle_track + +!=============================================================================== +! WRITE_PARTICLE_TRACK outputs particle position to a binary file. +! +! TODO: This subroutine needs to be modified to work with HDF5 files. Perhaps +! it should also be made somehow more general so that it can output +! information other than just particle position. +!=============================================================================== + + subroutine write_particle_track() + write(UNIT_TRACK) p % coord0 % xyz + end subroutine write_particle_track + +!=============================================================================== +! FINALIZE_PARTICLE_TRACK closes the particle track file. +! +! TODO: This subroutine needs to be modified to work with HDF5 files. +!=============================================================================== + + subroutine finalize_particle_track() + close(UNIT=UNIT_TRACK) + end subroutine finalize_particle_track + end module output diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 15afeb4cb5..39d8156fdb 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -11,7 +11,6 @@ module particle_restart use physics, only: transport use random_lcg, only: set_particle_seed use source, only: initialize_particle - use string, only: to_str #ifdef HDF5 use hdf5_interface @@ -20,7 +19,6 @@ module particle_restart implicit none private public :: run_particle_restart - public :: write_particle_track #ifdef HDF5 integer(HID_T) :: hdf5_particle_file @@ -144,36 +142,12 @@ contains current_gen - 1)*n_particles + p % id call set_particle_seed(particle_seed) - ! Open particle track output file. - if (write_track) then - filename = trim(path_output) // 'track_' // trim(to_str(current_batch)) & - // '_' // trim(to_str(current_work)) // '.binary' - open(UNIT=UNIT_TRACK, FILE=filename, STATUS='replace', & - ACCESS='stream') - end if - ! transport neutron call transport() - ! Close particle track output file. - if (write_track) then - close(UNIT=UNIT_TRACK) - endif - ! write output if particle made it call print_particle() end subroutine run_particle_restart -!=============================================================================== -! WRITE_PARTICLE_TRACK outputs particle position to a binary file. This -! subroutine needs to be modified to work with HDF5 files. Perhaps it should -! also be made somehow more general so that it can output information other -! than just particle position. -!=============================================================================== - - subroutine write_particle_track() - write(UNIT_TRACK) p % coord0 % xyz - end subroutine write_particle_track - end module particle_restart diff --git a/src/physics.F90 b/src/physics.F90 index ba8eac3fa5..a5f3318538 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -13,9 +13,11 @@ module physics use interpolation, only: interpolate_tab1 use material_header, only: Material use mesh, only: get_mesh_indices - use output, only: write_message + use output, only: write_message, initialize_particle_track, & + write_particle_track, & + finalize_particle_track use particle_header, only: LocalCoord - use particle_restart, only: write_particle_track +! use particle_restart, only: write_particle_track use particle_restart_write, only: write_particle_restart use random_lcg, only: prn use search, only: binary_search @@ -74,6 +76,11 @@ contains ! Force calculation of cross-sections by setting last energy to zero micro_xs % last_E = ZERO + ! Prepare to write out particle track. + if (write_track) then + call initialize_particle_track() + endif + do while (p % alive) ! Write particle track. @@ -106,9 +113,11 @@ contains coord % xyz = coord % xyz + distance * coord % uvw coord => coord % next end do + + ! Write particle track - if (run_mode == MODE_PARTICLE) then - write(UNIT_TRACK) p % coord0 % xyz + if (write_track) then + call write_particle_track() end if ! Score track-length tallies @@ -185,6 +194,11 @@ contains end do + ! Finish particle track. + if (write_track) then + call finalize_particle_track() + endif + end subroutine transport !=============================================================================== diff --git a/src/utils/track.py b/src/utils/track.py index a5b33b7258..1907d1ec0a 100755 --- a/src/utils/track.py +++ b/src/utils/track.py @@ -1,7 +1,19 @@ #!/usr/bin/env python2 """Convert binary particle track to VTK poly data. -Run 'track.py -help' for usage. +Usage information can be obtained by running 'track.py --help': + + usage: track.py [-h] [-o OUT] IN [IN ...] + + Convert binary particle track file to a .pvtp file. + + positional arguments: + IN Input particle track data filename(s). + + optional arguments: + -h, --help show this help message and exit + -o OUT, -out OUT, --out OUT + Output VTK poly data filename. """ @@ -17,13 +29,14 @@ def _parse_args(): description='Convert binary particle track file to a .pvtp file.') parser.add_argument('input', metavar='IN', type=argparse.FileType('rb'), nargs='+', - help='Input particle track data filename.') - parser.add_argument('-o', '-out', '--out', metavar='OUT', type=str, dest='out', + help='Input particle track data filename(s).') + parser.add_argument('-o', '-out', '--out', metavar='OUT', type=str, + dest='out', help='Output VTK poly data filename.') # Parse and return commandline arguments. return parser.parse_args() - + def main(): # Parse commandline argumetns. @@ -31,7 +44,7 @@ def main(): # Make sure that the output filename ends with '.pvtp'. if not args.out: - args.out = ''.join([os.path.splitext(args.input.name)[0], '.pvtp']) + args.out = 'tracks.pvtp' elif os.path.splitext(args.out)[1] != '.pvtp': args.out = ''.join([args.out, '.pvtp']) @@ -39,18 +52,14 @@ def main(): cells = vtk.vtkCellArray() j = 0 k = 0 - arr = vtk.vtkIntArray() for fin in args.input: track = fin.read() coords = [struct.unpack("ddd", track[24*i : 24*(i+1)]) for i in range(len(track)/24)] n_points = len(coords) - - i=0 + for triplet in coords: points.InsertNextPoint(triplet) - arr.InsertTuple1(i,k) - i+=1 # Create VTK line and assign points to line. line = vtk.vtkPolyLine() @@ -61,8 +70,7 @@ def main(): cells.InsertNextCell(line) j += n_points k += 1 - - points.SetData(arr) + global data data = vtk.vtkPolyData() data.SetPoints(points) data.SetLines(cells) @@ -78,9 +86,7 @@ if __name__ == '__main__': main() -#### The following code is retained for debugging purposes. If 'data' in the -#### in the function 'main' is made a global variable, the following code will -#### automatically plot the tracks. +#### The following code is retained for debugging purposes: ##poly_mapper = vtk.vtkPolyDataMapper() ##poly_mapper.SetInputConnection(data.GetProducerPort())