From 3571290a8f85d64d33f4fba30340fa91ed010996 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 11 Jan 2012 13:37:24 -0500 Subject: [PATCH 01/19] Improved error checking in cross_lattice subroutine. --- src/geometry.F90 | 47 ++++++++++++++++++++++++----------------------- 1 file changed, 24 insertions(+), 23 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 1dae6cedc8..70c0f8c306 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -495,29 +495,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 +506,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 From 01dc89d08199d7c561dfdb0ddef1899ad9c63b54 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 11 Jan 2012 13:43:19 -0500 Subject: [PATCH 02/19] Removed core_map utility. --- src/utils/core_map.py | 40 ---------------------------------------- 1 file changed, 40 deletions(-) delete mode 100755 src/utils/core_map.py 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() From 36173dc896e51ef945591008da40d488a2402346 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 12 Jan 2012 14:36:16 -0500 Subject: [PATCH 03/19] Added Andrew to list of advisors in documentation. --- docs/source/developers.rst | 2 ++ 1 file changed, 2 insertions(+) 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 From 069e795101377217a5f24a34ea4f1e2e2d94ac7f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 12 Jan 2012 15:57:35 -0500 Subject: [PATCH 04/19] Added option to turn probability tables off. --- src/input_xml.F90 | 3 +++ src/xml-fortran/templates/settings_t.xml | 1 + 2 files changed, 4 insertions(+) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e6c46cba60..db18cb3ab0 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -129,6 +129,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 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 @@ + From b207f271440ced5335ff12895a9033f216a179ab Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 17 Jan 2012 10:12:32 -0500 Subject: [PATCH 05/19] Many changes to Makefile. Made gnu default, added backtrace by default, added HDF5 linking information. --- src/Makefile | 92 +++++++++++++++++++++++++++++----------------------- 1 file changed, 52 insertions(+), 40 deletions(-) diff --git a/src/Makefile b/src/Makefile index 8b8a7b151e..c4f75e73ae 100644 --- a/src/Makefile +++ b/src/Makefile @@ -15,51 +15,25 @@ include OBJECTS # User Options #=============================================================================== -COMPILER = intel -DEBUG = no -PROFILE = no -OPTIMIZE = no -USE_MPI = no - -#=============================================================================== -# 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 = gnu +DEBUG = no +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 -Wall -fbacktrace LDFLAGS = # Debugging ifeq ($(DEBUG),yes) - F90FLAGS += -g -pedantic -std=f2008 -fbacktrace -fbounds-check \ + F90FLAGS += -g -pedantic -std=f2008 -fbounds-check \ -ffpe-trap=invalid,zero,overflow,underflow LDFLAGS += -g endif @@ -76,18 +50,45 @@ 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 + #=============================================================================== # PGI compiler options #=============================================================================== 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 @@ -151,11 +152,22 @@ endif # Use MPI for parallelism ifeq ($(USE_MPI),yes) - MPI = /opt/mpich2/1.4.1-intel + MPI_ROOT = /opt/mpich2/1.4.1-$(COMPILER) F90 = $(MPI)/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 #=============================================================================== @@ -166,7 +178,7 @@ ifeq ($(MACHINE),bluegene) endif #=============================================================================== -# Options for Cray XK6 ORNL Jaguar supercomputer +# Options for Cray XK6 ORNL Titan supercomputer #=============================================================================== ifeq ($(MACHINE),crayxk6) @@ -200,7 +212,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 From 6c0f36eedfbc3ec59f096212ce448ca99451afd6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 17 Jan 2012 12:44:30 -0500 Subject: [PATCH 06/19] Moved -Wall to DEBUG section for gfortran compiler in Makefile. --- src/Makefile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Makefile b/src/Makefile index c4f75e73ae..6efec7b337 100644 --- a/src/Makefile +++ b/src/Makefile @@ -28,12 +28,12 @@ USE_HDF5 = no ifeq ($(COMPILER),gnu) F90 = gfortran - F90FLAGS := -cpp -Wall -fbacktrace + F90FLAGS := -cpp -fbacktrace LDFLAGS = # Debugging ifeq ($(DEBUG),yes) - F90FLAGS += -g -pedantic -std=f2008 -fbounds-check \ + F90FLAGS += -g -Wall -pedantic -std=f2008 -fbounds-check \ -ffpe-trap=invalid,zero,overflow,underflow LDFLAGS += -g endif From 4bc81d501dc3ce1ba2c4200e502c305f9b007832 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 17 Jan 2012 13:13:39 -0500 Subject: [PATCH 07/19] Started HDF5 interface. --- src/DEPENDENCIES | 5 +++ src/OBJECTS | 1 + src/global.F90 | 11 +++++++ src/hdf5_interface.F90 | 72 ++++++++++++++++++++++++++++++++++++++++++ src/initialize.F90 | 21 +++++++++--- src/main.F90 | 11 +++++++ 6 files changed, 116 insertions(+), 5 deletions(-) create mode 100644 src/hdf5_interface.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 5016e4278b..bdcbc62d06 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -67,6 +67,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 @@ -77,6 +80,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 @@ -117,6 +121,7 @@ interpolation.o: string.o main.o: constants.o main.o: global.o +main.o: hdf5_interface.o main.o: initialize.o main.o: intercycle.o main.o: output.o diff --git a/src/OBJECTS b/src/OBJECTS index e7fc0db9bd..2cc565cf0e 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -14,6 +14,7 @@ fission.o \ geometry.o \ geometry_header.o \ global.o \ +hdf5_interface.o \ initialize.o \ intercycle.o \ interpolation.o \ diff --git a/src/global.F90 b/src/global.F90 index 5e8f5334d4..f7a5f482b2 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -17,6 +17,10 @@ module global use mpi #endif +#ifdef HDF5 + use hdf5 +#endif + implicit none save @@ -163,6 +167,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..d394c41973 --- /dev/null +++ b/src/hdf5_interface.F90 @@ -0,0 +1,72 @@ +module hdf5_interface + + use constants + use global + +#ifdef HDF5 + use hdf5, only: HSIZE_T, h5open_f, h5fcreate_f, h5fclose_f, h5close_f + use h5lt, only: h5ltmake_dataset_int_f +#endif + + implicit none + +contains + +#ifdef HDF5 + +!=============================================================================== +! OPEN_HDF5_OUTPUT +!=============================================================================== + + subroutine open_hdf5_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 open_hdf5_output + +!=============================================================================== +! WRITE_HDF5_SUMMARY +!=============================================================================== + + subroutine write_hdf5_summary() + + integer :: error ! Error flag + integer :: rank = 1 + integer(HSIZE_T) :: dims(1) = (/1/) + + ! 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) + + end subroutine write_hdf5_summary + +!=============================================================================== +! CLOSE_HDF5_OUTPUT +!=============================================================================== + + subroutine close_hdf5_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 close_hdf5_output + +#endif + +end module hdf5_interface diff --git a/src/initialize.F90 b/src/initialize.F90 index 1bcae2ef36..86e6b1dba4 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -25,6 +25,10 @@ module initialize use mpi #endif +#ifdef HDF5 + use hdf5_interface, only: open_hdf5_output +#endif + implicit none type(DictionaryII), pointer :: build_dict => null() @@ -50,13 +54,20 @@ 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 open_hdf5_output() +#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/main.F90 b/src/main.F90 index 199da60dc0..7d651ace14 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -19,6 +19,10 @@ program main use mpi #endif +#ifdef HDF5 + use hdf5_interface, only: write_hdf5_summary, close_hdf5_output +#endif + implicit none ! start timer for total run time @@ -45,6 +49,13 @@ program main ! deallocate arrays call free_memory() +#ifdef HDF5 + if (master) then + call write_hdf5_summary() + call close_hdf5_output() + end if +#endif + #ifdef MPI ! If MPI is in use and enabled, terminate it call MPI_FINALIZE(mpi_err) From e6aca9820a70aa83e54290e9b7a0aaec82366710 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 17 Jan 2012 16:21:13 -0500 Subject: [PATCH 08/19] Added finalize module to better group post-simulation subroutines. --- src/DEPENDENCIES | 7 ++++++- src/OBJECTS | 1 + src/finalize.F90 | 49 ++++++++++++++++++++++++++++++++++++++++++++++ src/initialize.F90 | 3 ++- src/main.F90 | 39 ++++++------------------------------ 5 files changed, 64 insertions(+), 35 deletions(-) create mode 100644 src/finalize.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index bdcbc62d06..d13f68d9dd 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -38,6 +38,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 @@ -121,7 +127,6 @@ interpolation.o: string.o main.o: constants.o main.o: global.o -main.o: hdf5_interface.o main.o: initialize.o main.o: intercycle.o main.o: output.o diff --git a/src/OBJECTS b/src/OBJECTS index 2cc565cf0e..f817f4cd5b 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -10,6 +10,7 @@ endf.o \ endf_header.o \ energy_grid.o \ error.o \ +finalize.o \ fission.o \ geometry.o \ geometry_header.o \ diff --git a/src/finalize.F90 b/src/finalize.F90 new file mode 100644 index 0000000000..7a4de9e86e --- /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: write_hdf5_summary, close_hdf5_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 write_hdf5_summary() + call close_hdf5_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/initialize.F90 b/src/initialize.F90 index 86e6b1dba4..c2361d599d 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -46,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 diff --git a/src/main.F90 b/src/main.F90 index 7d651ace14..c712746e6b 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -2,32 +2,25 @@ program main 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 use mpi #endif -#ifdef HDF5 - use hdf5_interface, only: write_hdf5_summary, close_hdf5_output -#endif - implicit none - ! start timer for total run time - call timer_start(time_total) - ! set up problem call initialize_run() @@ -36,30 +29,10 @@ program main 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() end if - - ! deallocate arrays - call free_memory() - -#ifdef HDF5 - if (master) then - call write_hdf5_summary() - call close_hdf5_output() - end if -#endif - -#ifdef MPI - ! If MPI is in use and enabled, terminate it - call MPI_FINALIZE(mpi_err) -#endif + + ! finalize run + call finalize_run() contains From e3bfbbe015d7788568153712f6cab9a2f03747db Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 17 Jan 2012 17:21:30 -0500 Subject: [PATCH 09/19] Fixed definition of F90 for MPI in Makefile. --- src/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Makefile b/src/Makefile index 6efec7b337..465845ff32 100644 --- a/src/Makefile +++ b/src/Makefile @@ -153,7 +153,7 @@ endif ifeq ($(USE_MPI),yes) MPI_ROOT = /opt/mpich2/1.4.1-$(COMPILER) - F90 = $(MPI)/bin/mpif90 + F90 = $(MPI_ROOT)/bin/mpif90 F90FLAGS += -DMPI endif From 3a4ac7996f784bc7cc92fe6796d57715bd9c1439 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 17 Jan 2012 17:22:12 -0500 Subject: [PATCH 10/19] Added date/time, MPI, and timing information in HDF5 output. --- src/finalize.F90 | 12 ++-- src/hdf5_interface.F90 | 123 +++++++++++++++++++++++++++++++++++------ src/initialize.F90 | 5 +- 3 files changed, 114 insertions(+), 26 deletions(-) diff --git a/src/finalize.F90 b/src/finalize.F90 index 7a4de9e86e..76ae8dc9a5 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -1,12 +1,12 @@ module finalize use global - use output, only: print_runtime - use tally, only: write_tallies, tally_statistics - use timing, only: timer_stop + use output, only: print_runtime + use tally, only: write_tallies, tally_statistics + use timing, only: timer_stop #ifdef HDF5 - use hdf5_interface, only: write_hdf5_summary, close_hdf5_output + use hdf5_interface, only: hdf5_write_timing, hdf5_close_output #endif implicit none @@ -31,8 +31,8 @@ contains #ifdef HDF5 if (master) then - call write_hdf5_summary() - call close_hdf5_output() + call hdf5_write_timing() + call hdf5_close_output() end if #endif diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index d394c41973..6ce0f01c40 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -4,8 +4,8 @@ module hdf5_interface use global #ifdef HDF5 - use hdf5, only: HSIZE_T, h5open_f, h5fcreate_f, h5fclose_f, h5close_f - use h5lt, only: h5ltmake_dataset_int_f + use hdf5 + use h5lt #endif implicit none @@ -15,10 +15,10 @@ contains #ifdef HDF5 !=============================================================================== -! OPEN_HDF5_OUTPUT +! HDF5_OPEN_OUTPUT !=============================================================================== - subroutine open_hdf5_output() + subroutine hdf5_open_output() character(9), parameter :: filename = "output.h5" ! File name integer :: error ! Error flag @@ -29,33 +29,120 @@ contains ! Create a new file using default properties. call h5fcreate_f(filename, H5F_ACC_TRUNC_F, hdf5_output_file, error) - end subroutine open_hdf5_output + end subroutine hdf5_open_output !=============================================================================== -! WRITE_HDF5_SUMMARY +! HDF5_WRITE_SUMMARY !=============================================================================== - subroutine write_hdf5_summary() + subroutine hdf5_write_summary() - integer :: error ! Error flag + 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) + 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) - end subroutine write_hdf5_summary + ! 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 !=============================================================================== -! CLOSE_HDF5_OUTPUT +! HDF5_WRITE_TIMING !=============================================================================== - subroutine close_hdf5_output() + 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 @@ -65,7 +152,7 @@ contains ! Close FORTRAN interface. call h5close_f(error) - end subroutine close_hdf5_output + end subroutine hdf5_close_output #endif diff --git a/src/initialize.F90 b/src/initialize.F90 index c2361d599d..d7791fc932 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -26,7 +26,7 @@ module initialize #endif #ifdef HDF5 - use hdf5_interface, only: open_hdf5_output + use hdf5_interface, only: hdf5_open_output, hdf5_write_summary #endif implicit none @@ -62,7 +62,8 @@ contains #ifdef HDF5 ! Open HDF5 output file for writing - call open_hdf5_output() + call hdf5_open_output() + call hdf5_write_summary() #endif ! Display title and initialization header From 97b8196c469e1598f7a74b777c24ce523a9b68b3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 18 Jan 2012 18:25:36 -0500 Subject: [PATCH 11/19] Small fix in print_particle subroutine. --- src/output.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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) From 7ba2b4d9acee431970005df3851a8d863a704bf9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 18 Jan 2012 18:26:16 -0500 Subject: [PATCH 12/19] Added check for zero cross section when calculating distance to collision. --- src/physics.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/physics.F90 b/src/physics.F90 index b0f1e2e2cd..68c701ae6a 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -79,7 +79,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) From c9ebf2424f34c34b6a5f3a8fb54ae3f7b86dfda3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 18 Jan 2012 18:34:01 -0500 Subject: [PATCH 13/19] Removed in_lower_universe attribute and changed test for reflection in lower universe. --- src/geometry.F90 | 4 +--- src/particle_header.F90 | 2 -- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 70c0f8c306..a006a72d93 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() 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 From 1571f275de6c923a07a267a42790f750e63a353c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 18 Jan 2012 18:38:28 -0500 Subject: [PATCH 14/19] Added precision check for distance to defined surface. --- src/geometry.F90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index a006a72d93..8fc16da83b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -851,10 +851,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 From 059784c2dbe9e62ea5ce258ba9df82a09676ea6a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 18 Jan 2012 19:26:31 -0500 Subject: [PATCH 15/19] Simplified handling of reflective boundary conditions. --- src/geometry.F90 | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 8fc16da83b..69e6747a2b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -305,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) @@ -323,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) @@ -336,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) @@ -349,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) @@ -362,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 @@ -378,14 +373,15 @@ 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 From 64a6c6d798971a30e6f7f10b78d625266e436897 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 18 Jan 2012 20:00:48 -0500 Subject: [PATCH 16/19] Improved distance to lattice calculation for cases with coincident surfaces. --- src/constants.F90 | 3 ++- src/geometry.F90 | 26 +++++++++++++------------- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index aae0fd6ba2..e8ff77d17c 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/geometry.F90 b/src/geometry.F90 index 69e6747a2b..ccbbb86dab 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -847,7 +847,7 @@ contains ! Check is calculated distance is new minimum if (d < dist) then - if (abs(d - dist)/dist >= FP_PRECISION) then + if (abs(d - dist)/dist >= FP_REL_PRECISION) then dist = d surface_crossed = -cl % surfaces(i) lattice_crossed = NONE @@ -869,16 +869,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, @@ -889,7 +889,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 @@ -901,16 +901,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 From f7d1fbb569cd20e54639ea548a5b02d00a59981c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 18 Jan 2012 20:19:55 -0500 Subject: [PATCH 17/19] Better handle reflection from surface coincident with lattice or universe boundary. --- src/geometry.F90 | 15 ++++++++++++++- src/physics.F90 | 1 + 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index ccbbb86dab..31b40045a2 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -386,6 +386,19 @@ contains 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 @@ -847,7 +860,7 @@ contains ! Check is calculated distance is new minimum if (d < dist) then - if (abs(d - dist)/dist >= FP_REL_PRECISION) then + if (abs(d - dist)/dist >= FP_PRECISION) then dist = d surface_crossed = -cl % surfaces(i) lattice_crossed = NONE diff --git a/src/physics.F90 b/src/physics.F90 index 68c701ae6a..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 From c6c2c07ff1827e83b3e6c5f5ca69ddb293b90e1c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Jan 2012 11:19:09 -0500 Subject: [PATCH 18/19] Check for invalid material density value. --- src/input_xml.F90 | 12 +++++++++--- src/xml-fortran/templates/materials_t.xml | 2 +- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index db18cb3ab0..c4a51e497a 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1,14 +1,14 @@ module input_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 @@ -496,7 +496,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)) 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 @@ - + From f3ae5d40ec30fb2e3fc6c4e28cd017c0dbd75377 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Jan 2012 11:27:56 -0500 Subject: [PATCH 19/19] Check for no directory in cross_sections.xml --- src/input_xml.F90 | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c4a51e497a..34ec0c6d25 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1032,8 +1032,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