diff --git a/docs/source/developers.rst b/docs/source/developers.rst index e26f5693eb..be48bb409b 100644 --- a/docs/source/developers.rst +++ b/docs/source/developers.rst @@ -13,8 +13,10 @@ Advisors to the project include: * `Benoit Forget`_ * `Kord Smith`_ +* `Andrew Siegel`_ .. _Paul Romano: mailto:romano7@mit.edu .. _Bryan Herman: mailto:bherman@mit.edu .. _Benoit Forget: mailto:bforget@mit.edu .. _Kord Smith: mailto:kord@mit.edu +.. _Andrew Siegel: mailto:siegela@mcs.anl.gov diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 0788a8d8de..1effb032d7 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -54,6 +54,12 @@ energy_grid.o: datatypes_header.o energy_grid.o: global.o energy_grid.o: output.o +finalize.o: global.o +finalize.o: hdf5_interface.o +finalize.o: output.o +finalize.o: tally.o +finalize.o: timing.o + fission.o: ace_header.o fission.o: constants.o fission.o: error.o @@ -84,6 +90,9 @@ global.o: source_header.o global.o: tally_header.o global.o: timing.o +hdf5_interface.o: constants.o +hdf5_interface.o: global.o + initialize.o: ace.o initialize.o: bank_header.o initialize.o: constants.o @@ -94,6 +103,7 @@ initialize.o: error.o initialize.o: geometry.o initialize.o: geometry_header.o initialize.o: global.o +initialize.o: hdf5_interface.o initialize.o: input_xml.o initialize.o: output.o initialize.o: random_lcg.o diff --git a/src/Makefile b/src/Makefile index a34e33b03c..ae362c5c7e 100644 --- a/src/Makefile +++ b/src/Makefile @@ -15,52 +15,25 @@ include OBJECTS # User Options #=============================================================================== -COMPILER = petsc -DEBUG = yes -PROFILE = no -OPTIMIZE = no -USE_MPI = no -USE_HDF5 = yes - -#=============================================================================== -# Intel Fortran compiler options -#=============================================================================== - -ifeq ($(COMPILER),intel) - F90 = ifort - F90FLAGS := -fpp -warn -assume byterecl - LDFLAGS = - - # Debugging - ifeq ($(DEBUG),yes) - F90FLAGS += -g -traceback -ftrapuv -fp-stack-check -check all -fpe0 - LDFLAGS += -g - endif - - # Profiling - ifeq ($(PROFILE),yes) - F90FLAGS += -pg - LDFLAGS += -pg - endif - - # Optimization - ifeq ($(OPTIMIZE),yes) - F90FLAGS += -O3 -ipo - endif -endif +COMPILER = yes +DEBUG = yes +PROFILE = no +OPTIMIZE = no +USE_MPI = no +USE_HDF5 = no #=============================================================================== # GNU Fortran compiler options #=============================================================================== -ifeq ($(COMPILER),gfortran) +ifeq ($(COMPILER),gnu) F90 = gfortran - F90FLAGS := -cpp -Wall + F90FLAGS := -cpp -fbacktrace LDFLAGS = # Debugging ifeq ($(DEBUG),yes) - F90FLAGS += -g -pedantic -std=f2008 -fbacktrace -fbounds-check \ + F90FLAGS += -g -Wall -pedantic -std=f2008 -fbounds-check \ -ffpe-trap=invalid,zero,overflow,underflow LDFLAGS += -g endif @@ -77,6 +50,33 @@ ifeq ($(COMPILER),gfortran) endif endif +#=============================================================================== +# Intel Fortran compiler options +#=============================================================================== + +ifeq ($(COMPILER),intel) + F90 = ifort + F90FLAGS := -cpp -warn -assume byterecl -traceback + LDFLAGS = + + # Debugging + ifeq ($(DEBUG),yes) + F90FLAGS += -g -ftrapuv -fp-stack-check -check all -fpe0 + LDFLAGS += -g + endif + + # Profiling + ifeq ($(PROFILE),yes) + F90FLAGS += -pg + LDFLAGS += -pg + endif + + # Optimization + ifeq ($(OPTIMIZE),yes) + F90FLAGS += -O3 -ipo + endif +endif + # use petsc compiler ifeq ($(COMPILER),petsc) F90 = ${FLINKER} @@ -88,11 +88,6 @@ ifeq ($(COMPILER),petsc) # LDFLAGS = ${PETSC_SYS_LIB} LDFLAGS = ${SLEPC_LIB} - ifeq ($(USE_HDF5),yes) - F90FLAGS += -I${HDF5_PREFIX}/include - LDFLAGS += -L${HDF5_PREFIX}/lib ${HDF5_PREFIX}/lib/libhdf5hl_fortran.a ${HDF5_PREFIX}/lib/libhdf5_hl.a ${HDF5_PREFIX}/lib/libhdf5_fortran.a ${HDF5_PREFIX}/lib/libhdf5.a -lz -lrt -lm -Wl,-rpath -Wl,${HDF5_PREFIX}/lib - endif - ifeq ($(DEBUG),yes) F90FLAGS += -g -fbacktrace -fbounds-check LDFLAGS += -g @@ -106,12 +101,12 @@ endif ifeq ($(COMPILER),pgi) F90 = pgf90 - F90FLAGS := -Mpreprocess -DNO_F2008 -Minform=inform + F90FLAGS := -Mpreprocess -DNO_F2008 -Minform=inform -traceback LDFLAGS = # Debugging ifeq ($(DEBUG),yes) - F90FLAGS += -g -Mbounds -Mchkptr -Mchkstk -traceback + F90FLAGS += -g -Mbounds -Mchkptr -Mchkstk LDFLAGS += -g endif @@ -175,11 +170,22 @@ endif # Use MPI for parallelism ifeq ($(USE_MPI),yes) - MPI = /opt/mpich2/1.4.1-intel - F90 = $(MPI)/bin/mpif90 + MPI_ROOT = /opt/mpich2/1.4.1-$(COMPILER) + F90 = $(MPI_ROOT)/bin/mpif90 F90FLAGS += -DMPI endif +# Use HDF5 for output + +ifeq ($(USE_HDF5),yes) + HDF5_ROOT = /opt/hdf5/1.8.8-$(COMPILER) + F90FLAGS += -DHDF5 -I${HDF5_ROOT}/include + LDFLAGS += -L${HDF5_ROOT}/lib ${HDF5_ROOT}/lib/libhdf5hl_fortran.a \ + ${HDF5_ROOT}/lib/libhdf5_hl.a ${HDF5_ROOT}/lib/libhdf5_fortran.a \ + ${HDF5_ROOT}/lib/libhdf5.a -lz -lrt -lm -Wl,-rpath -Wl,${HDF5_ROOT}/lib +endif + + #=============================================================================== # Options for IBM Blue Gene/P ANL supercomputer #=============================================================================== @@ -190,7 +196,7 @@ ifeq ($(MACHINE),bluegene) endif #=============================================================================== -# Options for Cray XK6 ORNL Jaguar supercomputer +# Options for Cray XK6 ORNL Titan supercomputer #=============================================================================== ifeq ($(MACHINE),crayxk6) @@ -224,7 +230,7 @@ neat: .PHONY: all xml-fortran clean neat distclean %.o: %.F90 - $(F90) -Ixml-fortran -Ixml-fortran/templates $(F90FLAGS) -c $< + $(F90) $(F90FLAGS) -Ixml-fortran -Ixml-fortran/templates -c $< #=============================================================================== # Dependencies diff --git a/src/OBJECTS b/src/OBJECTS index 43596c1b79..9f9f492c9d 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -13,10 +13,12 @@ endf.o \ endf_header.o \ energy_grid.o \ error.o \ +finalize.o \ fission.o \ geometry.o \ geometry_header.o \ global.o \ +hdf5_interface.o \ initialize.o \ intercycle.o \ interpolation.o \ diff --git a/src/constants.F90 b/src/constants.F90 index 20c59e79f2..58685a5677 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -22,7 +22,8 @@ module constants real(8), parameter :: TINY_BIT = 1e-8_8 ! User for precision in geometry - real(8), parameter :: FP_PRECISION = 1e-5_8 + real(8), parameter :: FP_PRECISION = 1e-14_8 + real(8), parameter :: FP_REL_PRECISION = 1e-5_8 ! Maximum number of collisions/crossings integer, parameter :: MAX_EVENTS = 10000 diff --git a/src/finalize.F90 b/src/finalize.F90 new file mode 100644 index 0000000000..76ae8dc9a5 --- /dev/null +++ b/src/finalize.F90 @@ -0,0 +1,49 @@ +module finalize + + use global + use output, only: print_runtime + use tally, only: write_tallies, tally_statistics + use timing, only: timer_stop + +#ifdef HDF5 + use hdf5_interface, only: hdf5_write_timing, hdf5_close_output +#endif + + implicit none + +contains + +!=============================================================================== +! FINALIZE_RUN +!=============================================================================== + + subroutine finalize_run() + + ! Calculate statistics for tallies and write to tallies.out + if (.not. plotting) then + call tally_statistics() + if (master) call write_tallies() + end if + + ! show timing statistics + call timer_stop(time_total) + if (master) call print_runtime() + +#ifdef HDF5 + if (master) then + call hdf5_write_timing() + call hdf5_close_output() + end if +#endif + + ! deallocate arrays + call free_memory() + +#ifdef MPI + ! If MPI is in use and enabled, terminate it + call MPI_FINALIZE(mpi_err) +#endif + + end subroutine finalize_run + +end module finalize diff --git a/src/geometry.F90 b/src/geometry.F90 index 1dae6cedc8..31b40045a2 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -149,7 +149,6 @@ contains ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL ! Create new level of coordinates - p % in_lower_universe = .true. allocate(p % coord % next) ! Move particle to next level and set universe @@ -187,7 +186,6 @@ contains else ! Create new level of coordinates - p % in_lower_universe = .true. allocate(p % coord % next) ! adjust local position of particle @@ -284,7 +282,7 @@ contains ! PARTICLE REFLECTS FROM SURFACE ! Do not handle reflective boundary conditions on lower universes - if (p % in_lower_universe) then + if (.not. associated(p % coord, p % coord0)) then message = "Cannot reflect particle " // trim(to_str(p % id)) // & " off surface in a lower universe." call fatal_error() @@ -307,11 +305,14 @@ contains select case (surf%type) case (SURF_PX) - p % coord0 % uvw = (/ -u, v, w /) + u = -u + case (SURF_PY) - p % coord0 % uvw = (/ u, -v, w /) + v = -v + case (SURF_PZ) - p % coord0 % uvw = (/ u, v, -w /) + w = -w + case (SURF_PLANE) ! Find surface coefficients and norm of vector normal to surface n1 = surf % coeffs(1) @@ -325,8 +326,6 @@ contains v = v - 2*dot_prod*n2/norm w = w - 2*dot_prod*n3/norm - ! Set vector - p % coord0 % uvw = (/ u, v, w /) case (SURF_CYL_X) ! Find y-y0, z-z0 and dot product of direction and surface normal y = p % coord0 % xyz(2) - surf % coeffs(1) @@ -338,8 +337,6 @@ contains v = v - 2*dot_prod*y/(R*R) w = w - 2*dot_prod*z/(R*R) - ! Set vector - p % coord0 % uvw = (/ u, v, w /) case (SURF_CYL_Y) ! Find x-x0, z-z0 and dot product of direction and surface normal x = p % coord0 % xyz(1) - surf % coeffs(1) @@ -351,8 +348,6 @@ contains u = u - 2*dot_prod*x/(R*R) w = w - 2*dot_prod*z/(R*R) - ! Set vector - p % coord0 % uvw = (/ u, v, w /) case (SURF_CYL_Z) ! Find x-x0, y-y0 and dot product of direction and surface normal x = p % coord0 % xyz(1) - surf % coeffs(1) @@ -364,8 +359,6 @@ contains u = u - 2*dot_prod*x/(R*R) v = v - 2*dot_prod*y/(R*R) - ! Set vector - p % coord0 % uvw = (/ u, v, w /) case (SURF_SPHERE) ! Find x-x0, y-y0, z-z0 and dot product of direction and surface ! normal @@ -380,18 +373,32 @@ contains v = v - 2*dot_prod*y/(R*R) w = w - 2*dot_prod*z/(R*R) - ! Set vector - p % coord0 % uvw = (/ u, v, w /) case default message = "Reflection not supported for surface " // & trim(to_str(surf % id)) call fatal_error() end select + ! Set new particle direction + p % coord0 % uvw = (/ u, v, w /) + ! Reassign particle's cell and surface p % coord0 % cell = last_cell p % surface = -p % surface + ! If a reflective surface is coincident with a lattice or universe + ! boundary, it is necessary to redetermine the particle's coordinates in + ! the lower universes. + + if (associated(p % coord0 % next)) then + call deallocate_coord(p % coord0 % next) + call find_cell(p, found) + if (.not. found) then + message = "Couldn't find particle after reflecting from surface." + call fatal_error() + end if + end if + ! Set previous coordinate going slightly past surface crossing p % last_xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw @@ -495,29 +502,7 @@ contains ! Check to make sure still in lattice i_x = p % coord % lattice_x i_y = p % coord % lattice_y - if (i_x < 1 .or. i_x > lat % n_x) then - message = "Particle " // trim(to_str(p % id)) // " reached edge of & - &lattice " // trim(to_str(lat % id)) // " at position (" // & - trim(to_str(i_x)) // "," // trim(to_str(i_y)) // ")." - call fatal_error() - elseif (i_y < 1 .or. i_y > lat % n_y) then - message = "Particle " // trim(to_str(p % id)) // " reached edge of & - &lattice " // trim(to_str(lat % id)) // " at position (" // & - trim(to_str(i_x)) // "," // trim(to_str(i_y)) // ")." - call fatal_error() - end if - - ! Find universe for next lattice element - p % coord % universe = lat % element(i_x, i_y) - - ! Find cell in next lattice element - call find_cell(p, found) - if (.not. found) then - ! In some circumstances, a particle crossing the corner of a cell may not - ! be able to be found in the next universe. In this scenario we cut off - ! all lower-level coordinates and search from universe zero - - ! Remove lower coordinates + if (i_x < 1 .or. i_x > lat % n_x .or. i_y < 1 .or. i_y > lat % n_y) then call deallocate_coord(p % coord0 % next) p % coord => p % coord0 @@ -528,6 +513,29 @@ contains " after crossing a lattice boundary." call fatal_error() end if + else + ! Find universe for next lattice element + p % coord % universe = lat % element(i_x, i_y) + + ! Find cell in next lattice element + call find_cell(p, found) + if (.not. found) then + ! In some circumstances, a particle crossing the corner of a cell may not + ! be able to be found in the next universe. In this scenario we cut off + ! all lower-level coordinates and search from universe zero + + ! Remove lower coordinates + call deallocate_coord(p % coord0 % next) + p % coord => p % coord0 + + ! Search for particle + call find_cell(p, found) + if (.not. found) then + message = "Could not locate particle " // trim(to_str(p % id)) // & + " after crossing a lattice boundary." + call fatal_error() + end if + end if end if end subroutine cross_lattice @@ -852,10 +860,12 @@ contains ! Check is calculated distance is new minimum if (d < dist) then - dist = d - surface_crossed = -cl % surfaces(i) - lattice_crossed = NONE - final_coord => coord + if (abs(d - dist)/dist >= FP_PRECISION) then + dist = d + surface_crossed = -cl % surfaces(i) + lattice_crossed = NONE + final_coord => coord + end if end if end do SURFACE_LOOP @@ -872,16 +882,16 @@ contains z = coord % xyz(3) ! determine oncoming edge - x0 = lat % width_x * 0.5_8 - y0 = lat % width_y * 0.5_8 + x0 = sign(lat % width_x * 0.5_8, u) + y0 = sign(lat % width_y * 0.5_8, v) ! left and right sides - if (u == ZERO) then + if (abs(x - x0) < FP_PRECISION) then + d = INFINITY + elseif (u == ZERO) then d = INFINITY - elseif (u > 0) then - d = (x0 - x)/u else - d = -(x + x0)/u + d = (x0 - x)/u end if ! If the lattice boundary is coincident with the parent cell boundary, @@ -892,7 +902,7 @@ contains ! point precision. if (d < dist) then - if (abs(d - dist)/dist >= FP_PRECISION) then + if (abs(d - dist)/dist >= FP_REL_PRECISION) then dist = d if (u > 0) then lattice_crossed = LATTICE_RIGHT @@ -904,16 +914,16 @@ contains end if ! top and bottom sides - if (v == ZERO) then + if (abs(y - y0) < FP_PRECISION) then + d = INFINITY + elseif (v == ZERO) then d = INFINITY - elseif (v > 0) then - d = (y0 - y)/v else - d = -(y + y0)/v + d = (y0 - y)/v end if if (d < dist) then - if (abs(d - dist)/dist >= FP_PRECISION) then + if (abs(d - dist)/dist >= FP_REL_PRECISION) then dist = d if (v > 0) then lattice_crossed = LATTICE_TOP diff --git a/src/global.F90 b/src/global.F90 index 93d0000488..44c10ef075 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -18,6 +18,10 @@ module global use mpi #endif +#ifdef HDF5 + use hdf5 +#endif + implicit none save @@ -164,6 +168,13 @@ module global real(8) :: plot_basis(6) real(8) :: pixel + ! ============================================================================ + ! HDF5 VARIABLES + +#ifdef HDF5 + integer(HID_T) :: hdf5_output_file +#endif + ! ============================================================================ ! MISCELLANEOUS VARIABLES diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 new file mode 100644 index 0000000000..6ce0f01c40 --- /dev/null +++ b/src/hdf5_interface.F90 @@ -0,0 +1,159 @@ +module hdf5_interface + + use constants + use global + +#ifdef HDF5 + use hdf5 + use h5lt +#endif + + implicit none + +contains + +#ifdef HDF5 + +!=============================================================================== +! HDF5_OPEN_OUTPUT +!=============================================================================== + + subroutine hdf5_open_output() + + character(9), parameter :: filename = "output.h5" ! File name + integer :: error ! Error flag + + ! Initialize FORTRAN interface. + call h5open_f (error) + + ! Create a new file using default properties. + call h5fcreate_f(filename, H5F_ACC_TRUNC_F, hdf5_output_file, error) + + end subroutine hdf5_open_output + +!=============================================================================== +! HDF5_WRITE_SUMMARY +!=============================================================================== + + subroutine hdf5_write_summary() + + integer :: error + integer :: rank = 1 + integer(HSIZE_T) :: dims(1) = (/1/) + character(8) :: date_ + character(10) :: time_ + character(19) :: current_time + + ! Write version information + call h5ltmake_dataset_int_f(hdf5_output_file, "version_major", & + rank, dims, (/ VERSION_MAJOR /), error) + call h5ltmake_dataset_int_f(hdf5_output_file, "version_minor", & + rank, dims, (/ VERSION_MINOR /), error) + call h5ltmake_dataset_int_f(hdf5_output_file, "version_release", & + rank, dims, (/ VERSION_RELEASE /), error) + + ! Write current date and time + call date_and_time(DATE=date_, TIME=time_) + current_time = date_(1:4) // "-" // date_(5:6) // "-" // date_(7:8) // & + " " // time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) + call h5ltmake_dataset_string_f(hdf5_output_file, "/date_and_time", & + current_time, error) + + ! Write MPI information + call h5ltmake_dataset_int_f(hdf5_output_file, "n_procs", & + rank, dims, (/ n_procs /), error) + call h5ltset_attribute_string_f(hdf5_output_file, "n_procs", & + "description", "Number of MPI processes", error) + + end subroutine hdf5_write_summary + +!=============================================================================== +! HDF5_WRITE_TIMING +!=============================================================================== + + subroutine hdf5_write_timing() + + integer :: error + integer :: rank = 1 + integer(HSIZE_T) :: dims(1) = (/1/) + integer(HID_T) :: timing_group + + ! Create group for timing + call h5gcreate_f(hdf5_output_file, "/timing", timing_group, error) + + ! Write timing data + call h5ltmake_dataset_double_f(timing_group, "time_initialize", & + rank, dims, (/ time_initialize % elapsed /), error) + call h5ltmake_dataset_double_f(timing_group, "time_read_xs", & + rank, dims, (/ time_read_xs % elapsed /), error) + call h5ltmake_dataset_double_f(timing_group, "time_unionize", & + rank, dims, (/ time_unionize % elapsed /), error) + call h5ltmake_dataset_double_f(timing_group, "time_compute", & + rank, dims, (/ time_compute % elapsed /), error) + call h5ltmake_dataset_double_f(timing_group, "time_intercycle", & + rank, dims, (/ time_intercycle % elapsed /), error) + call h5ltmake_dataset_double_f(timing_group, "time_tallies", & + rank, dims, (/ time_ic_tallies % elapsed /), error) + call h5ltmake_dataset_double_f(timing_group, "time_sample", & + rank, dims, (/ time_ic_sample % elapsed /), error) + call h5ltmake_dataset_double_f(timing_group, "time_sendrecv", & + rank, dims, (/ time_ic_sendrecv % elapsed /), error) + call h5ltmake_dataset_double_f(timing_group, "time_rebuild", & + rank, dims, (/ time_ic_rebuild % elapsed /), error) + call h5ltmake_dataset_double_f(timing_group, "time_inactive", & + rank, dims, (/ time_inactive % elapsed /), error) + call h5ltmake_dataset_double_f(timing_group, "time_active", & + rank, dims, (/ time_active % elapsed /), error) + call h5ltmake_dataset_double_f(timing_group, "time_total", & + rank, dims, (/ time_total % elapsed /), error) + + ! Add descriptions to timing data + call h5ltset_attribute_string_f(timing_group, "time_initialize", & + "description", "Total time elapsed for initialization (s)", error) + call h5ltset_attribute_string_f(timing_group, "time_read_xs", & + "description", "Time reading cross-section libraries (s)", error) + call h5ltset_attribute_string_f(timing_group, "time_unionize", & + "description", "Time unionizing energy grid (s)", error) + call h5ltset_attribute_string_f(timing_group, "time_compute", & + "description", "Total time in computation (s)", error) + call h5ltset_attribute_string_f(timing_group, "time_intercycle", & + "description", "Total time between cycles (s)", error) + call h5ltset_attribute_string_f(timing_group, "time_tallies", & + "description", "Time between cycles accumulating tallies (s)", error) + call h5ltset_attribute_string_f(timing_group, "time_sample", & + "description", "Time between cycles sampling source sites (s)", error) + call h5ltset_attribute_string_f(timing_group, "time_sendrecv", & + "description", "Time between cycles SEND/RECVing source sites (s)", error) + call h5ltset_attribute_string_f(timing_group, "time_rebuild", & + "description", "Time between cycles reconstructing source bank (s)", error) + call h5ltset_attribute_string_f(timing_group, "time_inactive", & + "description", "Total time in inactive cycles (s)", error) + call h5ltset_attribute_string_f(timing_group, "time_active", & + "description", "Total time in active cycles (s)", error) + call h5ltset_attribute_string_f(timing_group, "time_total", & + "description", "Total time elapsed (s)", error) + + ! Close timing group + call h5gclose_f(timing_group, error) + + end subroutine hdf5_write_timing + +!=============================================================================== +! HDF5_CLOSE_OUTPUT +!=============================================================================== + + subroutine hdf5_close_output() + + integer :: error ! Error flag + + ! Terminate access to the file. + call h5fclose_f(hdf5_output_file, error) + + ! Close FORTRAN interface. + call h5close_f(error) + + end subroutine hdf5_close_output + +#endif + +end module hdf5_interface diff --git a/src/initialize.F90 b/src/initialize.F90 index 1bcae2ef36..d7791fc932 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -25,6 +25,10 @@ module initialize use mpi #endif +#ifdef HDF5 + use hdf5_interface, only: hdf5_open_output, hdf5_write_summary +#endif + implicit none type(DictionaryII), pointer :: build_dict => null() @@ -42,7 +46,8 @@ contains type(Universe), pointer :: univ - ! Start initialization timer + ! Start total and initialization timer + call timer_start(time_total) call timer_start(time_initialize) ! Setup MPI @@ -50,13 +55,21 @@ contains ! Read command line arguments call read_command_line() - if (master) call create_summary_file() - ! Print the OpenMC title and version/date/time information - if (master) call title() + if (master) then + ! Create summary.out file + call create_summary_file() - ! Print initialization header block - if (master) call header("INITIALIZATION", level=1) +#ifdef HDF5 + ! Open HDF5 output file for writing + call hdf5_open_output() + call hdf5_write_summary() +#endif + + ! Display title and initialization header + call title() + call header("INITIALIZATION", level=1) + end if ! Initialize random number generator call initialize_prng() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index f395d92915..645ed3f899 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2,14 +2,14 @@ module input_xml use cmfd_utils, only: read_cmfd_xml use constants - use datatypes, only: dict_create, dict_add_key, dict_has_key, & + use datatypes, only: dict_create, dict_add_key, dict_has_key, & dict_get_key use error, only: fatal_error, warning use geometry_header, only: Cell, Surface, Lattice use global use mesh_header, only: StructuredMesh use output, only: write_message - use string, only: lower_case, to_str, str_to_int, str_to_real, & + use string, only: lower_case, to_str, str_to_int, str_to_real, & split_string, starts_with, ends_with use tally_header, only: TallyObject @@ -131,6 +131,9 @@ contains ! Survival biasing if (trim(survival_) == 'on') survival_biasing = .true. + ! Probability tables + if (ptables_ == 'off') urr_ptables_on = .false. + ! Cutoffs if (size(cutoff_) > 0) then weight_cutoff = cutoff_(1) % weight @@ -495,7 +498,13 @@ contains ! Copy density -- the default value for the units is given in the ! material_t.xml file and doesn't need to be specified here, hence case ! default results in an error. - val = material_(i) % density % value + val = material_(i) % density % value + if (val <= ZERO) then + message = "Need to specify a positive density on material " // & + trim(to_str(m % id)) // "." + call fatal_error() + end if + units = material_(i) % density % units call lower_case(units) select case(trim(units)) @@ -1025,8 +1034,15 @@ contains ! Parse cross_sections.xml file call read_xml_file_cross_sections_t(path_cross_sections) - ! Copy directory information if present - directory = trim(directory_) + if (len_trim(directory_) > 0) then + ! Copy directory information if present + directory = trim(directory_) + else + ! If no directory is listed in cross_sections.xml, by default select the + ! directory in which the cross_sections.xml file resides + i = index(path_cross_sections, "/", BACK=.true.) + directory = path_cross_sections(1:i) + end if ! determine whether binary/ascii if (filetype_ == 'ascii') then diff --git a/src/main.F90 b/src/main.F90 index 63a89c3f9a..1ec117f2e5 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -4,17 +4,17 @@ program main use cmfd_utils, only: read_hdf5 use constants use global + use finalize, only: finalize_run use initialize, only: initialize_run use intercycle, only: shannon_entropy, calculate_keff, synchronize_bank - use output, only: write_message, header, print_runtime + use output, only: write_message, header use particle_header, only: Particle use plot, only: run_plot use physics, only: transport use random_lcg, only: set_particle_seed use source, only: get_source_particle use string, only: to_str - use tally, only: synchronize_tallies, write_tallies, & - tally_statistics + use tally, only: synchronize_tallies use timing, only: timer_start, timer_stop #ifdef MPI @@ -23,6 +23,7 @@ program main implicit none +<<<<<<< HEAD integer :: nargs ! number of arguments character(len=32) :: arg ! the argument @@ -57,32 +58,18 @@ program main case default - ! start timer for total run time - call timer_start(time_total) - ! set up problem call initialize_run() ! start problem if (plotting) then - call run_plot() + call run_plot() else - call run_problem() - - ! Calculate statistics for tallies and write to tallies.out - call tally_statistics() - if (master) call write_tallies() - - ! show timing statistics - call timer_stop(time_total) - if (master) call print_runtime() + call run_problem() end if - ! call cmfd calculation - call execute_cmfd() - - ! deallocate arrays - call free_memory() + ! finalize run + call finalize_run() end select diff --git a/src/output.F90 b/src/output.F90 index 7b2898528c..4acbf169a2 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -261,7 +261,7 @@ contains ! Print surface if (p % surface /= NONE) then s => surfaces(p % surface) - write(ou,*) ' Surface = ' // to_str(s % id) + write(ou,*) ' Surface = ' // to_str(s % id) end if write(ou,*) ' Weight = ' // to_str(p % wgt) diff --git a/src/particle_header.F90 b/src/particle_header.F90 index fcd6068039..8f81fbeb1f 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -39,7 +39,6 @@ module particle_header integer :: type ! Particle type (n, p, e, etc) ! Particle coordinates - logical :: in_lower_universe ! is particle in lower universe? type(LocalCoord), pointer :: coord0 => null() ! coordinates on universe 0 type(LocalCoord), pointer :: coord => null() ! coordinates on lowest universe @@ -107,7 +106,6 @@ contains allocate(p % coord0) p % coord0 % universe = BASE_UNIVERSE p % coord => p % coord0 - p % in_lower_universe = .false. end subroutine initialize_particle diff --git a/src/physics.F90 b/src/physics.F90 index b0f1e2e2cd..2d07d2f4e9 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -45,6 +45,7 @@ contains ! Particle couldn't be located if (.not. found_cell) then message = "Could not locate particle " // trim(to_str(p % id)) + print *, p % coord0 % xyz call fatal_error() end if @@ -79,7 +80,11 @@ contains call distance_to_boundary(p, d_boundary, surface_crossed, lattice_crossed) ! Sample a distance to collision - d_collision = -log(prn()) / material_xs % total + if (material_xs % total == ZERO) then + d_collision = INFINITY + else + d_collision = -log(prn()) / material_xs % total + end if ! Select smaller of the two distances distance = min(d_boundary, d_collision) diff --git a/src/utils/core_map.py b/src/utils/core_map.py deleted file mode 100755 index c694f79775..0000000000 --- a/src/utils/core_map.py +++ /dev/null @@ -1,40 +0,0 @@ -#!/usr/bin/env python - -import sys - -import numpy as np -import matplotlib.pyplot as pyplot - -# Get filename and file object -filename = sys.argv[1] -axial_level = int(sys.argv[2]) -if len(sys.argv) > 3: - stage = int(sys.argv[3]) -else: - stage = 1 -fh = open(filename,'r') - -# Read size of mesh -words = fh.readline().split() -nx, ny, nz, nstage = [int(item) for item in words] - -# Read values -value_dict = {} -for line in fh: - words = line.split() - i, j, k, m = [int(item) for item in words[:4]] - value = float(words[-1]) - value_dict[i,j,k,m] = value - -# Set up matrix -matrix = np.array([[value_dict[i+1,j+1,axial_level,stage] for i in range(nx)] - for j in range(ny)]) - -# Make figure -pyplot.pcolor(matrix) -pyplot.colorbar() -pyplot.xlim([0,nx]) -pyplot.ylim([0,ny]) -pyplot.xticks([]) -pyplot.yticks([]) -pyplot.show() diff --git a/src/xml-fortran/templates/materials_t.xml b/src/xml-fortran/templates/materials_t.xml index 5e5a0ee41f..45225412e2 100644 --- a/src/xml-fortran/templates/materials_t.xml +++ b/src/xml-fortran/templates/materials_t.xml @@ -8,7 +8,7 @@ - + diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml index c69439612e..0bfa90278c 100644 --- a/src/xml-fortran/templates/settings_t.xml +++ b/src/xml-fortran/templates/settings_t.xml @@ -23,6 +23,7 @@ +