From 6987cb18071598d7a0456eab097356b7497f3bf6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 11 Dec 2011 15:22:46 -0500 Subject: [PATCH 01/33] Moved summary output information into summary.out. Removed logging module. --- src/DEPENDENCIES | 4 - src/OBJECTS | 1 - src/constants.F90 | 6 +- src/initialize.F90 | 7 +- src/logging.F90 | 43 ------ src/main.F90 | 4 +- src/output.F90 | 371 +++++++++++++++++++++++++++------------------ 7 files changed, 229 insertions(+), 207 deletions(-) delete mode 100644 src/logging.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index c1a0c2c08f..f6585ba4dc 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -72,7 +72,6 @@ initialize.o: geometry.o initialize.o: geometry_header.o initialize.o: global.o initialize.o: input_xml.o -initialize.o: logging.o initialize.o: mpi_routines.o initialize.o: output.o initialize.o: random_lcg.o @@ -103,9 +102,6 @@ interpolation.o: global.o interpolation.o: search.o interpolation.o: string.o -logging.o: constants.o -logging.o: global.o - main.o: constants.o main.o: global.o main.o: initialize.o diff --git a/src/OBJECTS b/src/OBJECTS index 6ce0dcafc7..e8df9f0f7e 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -17,7 +17,6 @@ global.o \ initialize.o \ interpolation.o \ input_xml.o \ -logging.o \ main.o \ material_header.o \ mesh_header.o \ diff --git a/src/constants.F90 b/src/constants.F90 index eba9786d66..352d500c4b 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -280,9 +280,9 @@ module constants PROB_CRITICALITY = 2 ! Criticality problem ! Unit numbers - integer, parameter :: UNIT_LOG = 11 ! unit # for writing log file - integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file - integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot 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 end module constants diff --git a/src/initialize.F90 b/src/initialize.F90 index cc5314d125..621df2c649 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -12,10 +12,9 @@ module initialize use global use input_xml, only: read_input_xml, read_cross_sections_xml, & cells_in_univ_dict - use logging, only: create_log use mpi_routines, only: setup_mpi use output, only: title, header, print_summary, print_geometry, & - print_plot + print_plot, create_summary_file use random_lcg, only: initialize_prng use source, only: initialize_source use string, only: int_to_str, starts_with, ends_with, lower_case @@ -47,13 +46,13 @@ contains ! Read command line arguments call read_command_line() - if (master) call create_log() + if (master) call create_summary_file() ! Print the OpenMC title and version/date/time information if (master) call title() ! Print initialization header block - if (master) call header("INITIALIZATION", 1) + if (master) call header("INITIALIZATION", level=1) ! Initialize random number generator call initialize_prng() diff --git a/src/logging.F90 b/src/logging.F90 deleted file mode 100644 index 377c1ab827..0000000000 --- a/src/logging.F90 +++ /dev/null @@ -1,43 +0,0 @@ -module logging - - use constants, only: MAX_WORD_LEN, UNIT_LOG - use global, only: path_input - - implicit none - -contains - -!=============================================================================== -! CREATE_LOG creates a new log file (or overwrites the existing log) -!=============================================================================== - - subroutine create_log() - - character(MAX_WORD_LEN) :: path_log ! path of log file - logical :: file_exists ! does log file already exist? - ! integer :: ioError ! error status for file access - - ! Create filename for log file - path_log = trim(path_input) // ".log" - - ! Check if log file already exists - inquire(FILE=path_log, EXIST=file_exists) - if (file_exists) then - ! Possibly copy old log file - end if - - ! Open log file for writing - ! open(FILE=path_log, UNIT=UNIT_LOG, STATUS='replace', & - ! & ACTION='write', IOSTAT=ioError) - - end subroutine create_log - -!=============================================================================== -! LOG_TALLIES -!=============================================================================== - - subroutine log_tallies() - - end subroutine log_tallies - -end module logging diff --git a/src/main.F90 b/src/main.F90 index a123d4790c..35479e9ea1 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -58,7 +58,7 @@ contains integer(8) :: i_particle ! history index type(Particle), pointer :: p => null() - if (master) call header("BEGIN SIMULATION", 1) + if (master) call header("BEGIN SIMULATION", level=1) tallies_on = .false. call timer_start(time_inactive) @@ -142,7 +142,7 @@ contains ! ========================================================================== ! END OF RUN WRAPUP - if (master) call header("SIMULATION FINISHED", 1) + if (master) call header("SIMULATION FINISHED", level=1) end subroutine run_problem diff --git a/src/output.F90 b/src/output.F90 index f0f7d6c9a5..371f239474 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -31,31 +31,38 @@ contains character(10) :: today_date character(8) :: today_time - write(ou,*) - write(ou,*) ' .d88888b. 888b d888 .d8888b.' - write(ou,*) ' d88P" "Y88b 8888b d8888 d88P Y88b' - write(ou,*) ' 888 888 88888b.d88888 888 888' - write(ou,*) ' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ' - write(ou,*) ' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ' - write(ou,*) ' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888' - write(ou,*) ' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P' - write(ou,*) ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"' - write(ou,*) '__________________888______________________________________________________' - write(ou,*) ' 888' - write(ou,*) ' 888' - write(ou,*) + write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') & + ' .d88888b. 888b d888 .d8888b.', & + ' d88P" "Y88b 8888b d8888 d88P Y88b', & + ' 888 888 88888b.d88888 888 888', & + ' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', & + ' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', & + ' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', & + ' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', & + ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', & + '__________________888______________________________________________________', & + ' 888', & + ' 888' ! Write version information - write(ou,*) ' Developed At: Massachusetts Institute of Technology' - write(ou,*) ' Lead Developer: Paul K. Romano' - write(ou,100) VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE -100 format (6X,"Version:",9X,I1,".",I1,".",I1) + write(UNIT=OUTPUT_UNIT, FMT=*) & + ' Developed At: Massachusetts Institute of Technology' + write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') & + VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE ! Write the date and time call get_today(today_date, today_time) - write(ou,101) trim(today_date), trim(today_time) -101 format (6X,"Date/Time:",7X,A,1X,A) - write(ou,*) + write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A,1X,A/)') & + trim(today_date), trim(today_time) + + ! Write information to summary file + call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1) + write(UNIT=UNIT_SUMMARY, FMT=*) & + "Copyright: 2010 Massachusetts Institute of Technology" + write(UNIT=UNIT_SUMMARY, FMT='(1X,A,5X,2(I1,"."),I1)') & + "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE + write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",3X,A,1X,A)') & + trim(today_date), trim(today_time) end subroutine title @@ -64,34 +71,32 @@ contains ! specified, it is assumed to be a minor header block (H3). !=============================================================================== - subroutine header(msg, level, unit_file) + subroutine header(msg, unit, level) character(*), intent(in) :: msg - integer, optional :: level - integer, optional :: unit_file + integer, optional :: unit + integer, optional :: level + integer :: n + integer :: m + integer :: unit_ integer :: header_level - integer :: n, m - integer :: unit - character(75) :: line + character(MAX_LINE_LEN) :: line - ! set default header level - if (.not. present(level)) then - header_level = 3 - else + ! set default level + if (present(level)) then header_level = level + else + header_level = 3 end if ! set default unit - if (present(unit_file)) then - unit = unit_file + if (present(unit)) then + unit_ = unit else - unit = ou + unit_ = OUTPUT_UNIT end if - ! Print first blank line - write(unit,*) - ! determine how many times to repeat '=' character n = (63 - len_trim(msg))/2 m = n @@ -104,23 +109,16 @@ contains ! print header based on level select case (header_level) case (1) - ! determine number of spaces to put in from of header - write(unit,*) repeat('=', 75) - write(unit,*) repeat('=', n) // '> ' // trim(line) // ' <' // & - & repeat('=', m) - write(unit,*) repeat('=', 75) + write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), & + repeat('=', n) // '> ' // trim(line) // ' <' // & + repeat('=', m), repeat('=', 75) case (2) - write(unit,*) trim(line) - write(unit,*) repeat('-', 75) + write(UNIT=unit_, FMT='(/2(1X,A/))') trim(line), repeat('-', 75) case (3) - n = (63 - len_trim(line))/2 - write(unit,*) repeat('=', n) // '> ' // trim(line) // ' <' // & - & repeat('=', m) + write(UNIT=unit_, FMT='(/1X,A/)') repeat('=', n) // '> ' // & + trim(line) // ' <' // repeat('=', m) end select - ! Print trailing blank line - write(unit, *) - end subroutine header !=============================================================================== @@ -290,39 +288,47 @@ contains ! PRINT_CELL displays the attributes of a cell !=============================================================================== - subroutine print_cell(c) + subroutine print_cell(c, unit) type(Cell), pointer :: c + integer, optional :: unit - integer :: temp - integer :: i + integer :: temp + integer :: i + integer :: unit_ character(MAX_LINE_LEN) :: string type(Universe), pointer :: u => null() type(Lattice), pointer :: l => null() type(Material), pointer :: m => null() - write(ou,*) 'Cell ' // int_to_str(c % id) + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + write(unit_,*) 'Cell ' // int_to_str(c % id) temp = dict_get_key(cell_dict, c % id) - write(ou,*) ' Array Index = ' // int_to_str(temp) + write(unit_,*) ' Array Index = ' // int_to_str(temp) u => universes(c % universe) - write(ou,*) ' Universe = ' // int_to_str(u % id) + write(unit_,*) ' Universe = ' // int_to_str(u % id) select case (c % type) case (CELL_NORMAL) - write(ou,*) ' Fill = NONE' + write(unit_,*) ' Fill = NONE' case (CELL_FILL) u => universes(c % fill) - write(ou,*) ' Fill = Universe ' // int_to_str(u % id) + write(unit_,*) ' Fill = Universe ' // int_to_str(u % id) case (CELL_LATTICE) l => lattices(c % fill) - write(ou,*) ' Fill = Lattice ' // int_to_str(l % id) + write(unit_,*) ' Fill = Lattice ' // int_to_str(l % id) end select if (c % material == 0) then - write(ou,*) ' Material = NONE' + write(unit_,*) ' Material = NONE' else m => materials(c % material) - write(ou,*) ' Material = ' // int_to_str(m % id) + write(unit_,*) ' Material = ' // int_to_str(m % id) end if - write(ou,*) ' Parent Cell = ' // int_to_str(c % parent) + write(unit_,*) ' Parent Cell = ' // int_to_str(c % parent) string = "" do i = 1, c % n_surfaces select case (c % surfaces(i)) @@ -338,8 +344,8 @@ contains string = trim(string) // ' ' // int_to_str(c % surfaces(i)) end select end do - write(ou,*) ' Surface Specification:' // trim(string) - write(ou,*) + write(unit_,*) ' Surface Specification:' // trim(string) + write(unit_,*) end subroutine print_cell @@ -347,23 +353,31 @@ contains ! PRINT_UNIVERSE displays the attributes of a universe !=============================================================================== - subroutine print_universe(univ) + subroutine print_universe(univ, unit) type(Universe), pointer :: univ + integer, optional :: unit - integer :: i + integer :: i + integer :: unit_ character(MAX_LINE_LEN) :: string type(Cell), pointer :: c => null() - write(ou,*) 'Universe ' // int_to_str(univ % id) - write(ou,*) ' Level = ' // int_to_str(univ % level) + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + write(unit_,*) 'Universe ' // int_to_str(univ % id) + write(unit_,*) ' Level = ' // int_to_str(univ % level) string = "" do i = 1, univ % n_cells c => cells(univ % cells(i)) string = trim(string) // ' ' // int_to_str(c % id) end do - write(ou,*) ' Cells =' // trim(string) - write(ou,*) + write(unit_,*) ' Cells =' // trim(string) + write(unit_,*) end subroutine print_universe @@ -371,18 +385,27 @@ contains ! PRINT_LATTICE displays the attributes of a lattice !=============================================================================== - subroutine print_lattice(lat) + subroutine print_lattice(lat, unit) type(Lattice), pointer :: lat + integer, optional :: unit - write(ou,*) 'Lattice ' // int_to_str(lat % id) - write(ou,*) ' n_x = ' // int_to_str(lat % n_x) - write(ou,*) ' n_y = ' // int_to_str(lat % n_y) - write(ou,*) ' x0 = ' // real_to_str(lat % x0) - write(ou,*) ' y0 = ' // real_to_str(lat % y0) - write(ou,*) ' width_x = ' // real_to_str(lat % width_x) - write(ou,*) ' width_y = ' // real_to_str(lat % width_y) - write(ou,*) + integer :: unit_ + + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + write(unit_,*) 'Lattice ' // int_to_str(lat % id) + write(unit_,*) ' n_x = ' // int_to_str(lat % n_x) + write(unit_,*) ' n_y = ' // int_to_str(lat % n_y) + write(unit_,*) ' x0 = ' // real_to_str(lat % x0) + write(unit_,*) ' y0 = ' // real_to_str(lat % y0) + write(unit_,*) ' width_x = ' // real_to_str(lat % width_x) + write(unit_,*) ' width_y = ' // real_to_str(lat % width_y) + write(unit_,*) end subroutine print_lattice @@ -390,14 +413,22 @@ contains ! PRINT_SURFACE displays the attributes of a surface !=============================================================================== - subroutine print_surface(surf) + subroutine print_surface(surf, unit) type(Surface), pointer :: surf + integer, optional :: unit integer :: i + integer :: unit_ character(MAX_LINE_LEN) :: string - write(ou,*) 'Surface ' // int_to_str(surf % id) + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + write(unit_,*) 'Surface ' // int_to_str(surf % id) select case (surf % type) case (SURF_PX) string = "X Plane" @@ -422,13 +453,13 @@ contains case (SURF_GQ) string = "General Quadratic" end select - write(ou,*) ' Type = ' // trim(string) + write(unit_,*) ' Type = ' // trim(string) string = "" do i = 1, size(surf % coeffs) string = trim(string) // ' ' // real_to_str(surf % coeffs(i), 4) end do - write(ou,*) ' Coefficients = ' // trim(string) + write(unit_,*) ' Coefficients = ' // trim(string) string = "" if (allocated(surf % neighbor_pos)) then @@ -436,7 +467,7 @@ contains string = trim(string) // ' ' // int_to_str(surf % neighbor_pos(i)) end do end if - write(ou,*) ' Positive Neighbors = ' // trim(string) + write(unit_,*) ' Positive Neighbors = ' // trim(string) string = "" if (allocated(surf % neighbor_neg)) then @@ -444,18 +475,18 @@ contains string = trim(string) // ' ' // int_to_str(surf % neighbor_neg(i)) end do end if - write(ou,*) ' Negative Neighbors =' // trim(string) + write(unit_,*) ' Negative Neighbors =' // trim(string) select case (surf % bc) case (BC_TRANSMIT) - write(ou,*) ' Boundary Condition = Transmission' + write(unit_,*) ' Boundary Condition = Transmission' case (BC_VACUUM) - write(ou,*) ' Boundary Condition = Vacuum' + write(unit_,*) ' Boundary Condition = Vacuum' case (BC_REFLECT) - write(ou,*) ' Boundary Condition = Reflective' + write(unit_,*) ' Boundary Condition = Reflective' case (BC_PERIODIC) - write(ou,*) ' Boundary Condition = Periodic' + write(unit_,*) ' Boundary Condition = Periodic' end select - write(ou,*) + write(unit_,*) end subroutine print_surface @@ -463,37 +494,45 @@ contains ! PRINT_MATERIAL displays the attributes of a material !=============================================================================== - subroutine print_material(mat) + subroutine print_material(mat, unit) type(Material), pointer :: mat + integer, optional :: unit - integer :: i - real(8) :: density + integer :: i + integer :: unit_ + real(8) :: density character(MAX_LINE_LEN) :: string type(Nuclide), pointer :: nuc => null() + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + ! Write identifier for material - write(ou,*) 'Material ' // int_to_str(mat % id) + write(unit_,*) 'Material ' // int_to_str(mat % id) ! Write total atom density in atom/b-cm - write(ou,*) ' Atom Density = ' // trim(real_to_str(mat % density)) & + write(unit_,*) ' Atom Density = ' // trim(real_to_str(mat % density)) & & // ' atom/b-cm' ! Write atom density for each nuclide in material - write(ou,*) ' Nuclides:' + write(unit_,*) ' Nuclides:' do i = 1, mat % n_nuclides nuc => nuclides(mat % nuclide(i)) density = mat % atom_density(i) string = ' ' // trim(nuc % name) // ' = ' // & & trim(real_to_str(density)) // ' atom/b-cm' - write(ou,*) trim(string) + write(unit_,*) trim(string) end do ! Write information on S(a,b) table if (mat % has_sab_table) then - write(ou,*) ' S(a,b) table = ' // trim(mat % sab_name) + write(unit_,*) ' S(a,b) table = ' // trim(mat % sab_name) end if - write(ou,*) + write(unit_,*) end subroutine print_material @@ -501,20 +540,28 @@ contains ! PRINT_TALLY displays the attributes of a tally !=============================================================================== - subroutine print_tally(t) + subroutine print_tally(t, unit) type(TallyObject), pointer :: t + integer, optional :: unit - integer :: i - integer :: id - character(MAX_LINE_LEN) :: string + integer :: i + integer :: id + integer :: unit_ + character(MAX_LINE_LEN) :: string type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() type(Universe), pointer :: u => null() type(Material), pointer :: m => null() type(StructuredMesh), pointer :: sm => null() - write(ou,*) 'Tally ' // int_to_str(t % id) + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + write(unit_,*) 'Tally ' // int_to_str(t % id) if (t % n_bins(T_CELL) > 0) then string = "" @@ -523,7 +570,7 @@ contains c => cells(id) string = trim(string) // ' ' // trim(int_to_str(c % id)) end do - write(ou, *) ' Cell Bins:' // trim(string) + write(unit_, *) ' Cell Bins:' // trim(string) end if if (t % n_bins(T_SURFACE) > 0) then @@ -533,7 +580,7 @@ contains s => surfaces(id) string = trim(string) // ' ' // trim(int_to_str(s % id)) end do - write(ou, *) ' Surface Bins:' // trim(string) + write(unit_, *) ' Surface Bins:' // trim(string) end if if (t % n_bins(T_UNIVERSE) > 0) then @@ -543,7 +590,7 @@ contains u => universes(id) string = trim(string) // ' ' // trim(int_to_str(u % id)) end do - write(ou, *) ' Material Bins:' // trim(string) + write(unit_, *) ' Material Bins:' // trim(string) end if if (t % n_bins(T_MATERIAL) > 0) then @@ -553,7 +600,7 @@ contains m => materials(id) string = trim(string) // ' ' // trim(int_to_str(m % id)) end do - write(ou, *) ' Material Bins:' // trim(string) + write(unit_, *) ' Material Bins:' // trim(string) end if if (t % n_bins(T_MESH) > 0) then @@ -564,7 +611,7 @@ contains do i = 2, sm % n_dimension string = trim(string) // ' x ' // trim(int_to_str(sm % dimension(i))) end do - write(ou, *) ' Mesh Bins:' // trim(string) + write(unit_, *) ' Mesh Bins:' // trim(string) end if if (t % n_bins(T_CELLBORN) > 0) then @@ -574,7 +621,7 @@ contains c => cells(id) string = trim(string) // ' ' // trim(int_to_str(c % id)) end do - write(ou, *) ' Birth Region Bins:' // trim(string) + write(unit_, *) ' Birth Region Bins:' // trim(string) end if if (t % n_bins(T_ENERGYIN) > 0) then @@ -583,7 +630,7 @@ contains string = trim(string) // ' ' // trim(real_to_str(& t % energy_in(i))) end do - write(ou,*) ' Incoming Energy Bins:' // trim(string) + write(unit_,*) ' Incoming Energy Bins:' // trim(string) end if if (t % n_bins(T_ENERGYOUT) > 0) then @@ -592,7 +639,7 @@ contains string = trim(string) // ' ' // trim(real_to_str(& t % energy_out(i))) end do - write(ou,*) ' Outgoing Energy Bins:' // trim(string) + write(unit_,*) ' Outgoing Energy Bins:' // trim(string) end if if (t % n_macro_bins > 0) then @@ -613,9 +660,9 @@ contains string = trim(string) // ' nu-fission' end select end do - write(ou,*) ' Macro Reactions:' // trim(string) + write(unit_,*) ' Macro Reactions:' // trim(string) end if - write(ou,*) + write(unit_,*) end subroutine print_tally @@ -632,36 +679,36 @@ contains type(Universe), pointer :: u => null() type(Lattice), pointer :: l => null() + ! print summary of surfaces + call header("SURFACE SUMMARY", unit=UNIT_SUMMARY) + do i = 1, n_surfaces + s => surfaces(i) + call print_surface(s, unit=UNIT_SUMMARY) + end do + ! print summary of cells - call header("CELL SUMMARY") + call header("CELL SUMMARY", unit=UNIT_SUMMARY) do i = 1, n_cells c => cells(i) - call print_cell(c) + call print_cell(c, unit=UNIT_SUMMARY) end do ! print summary of universes - call header("UNIVERSE SUMMARY") + call header("UNIVERSE SUMMARY", unit=UNIT_SUMMARY) do i = 1, n_universes u => universes(i) - call print_universe(u) + call print_universe(u, unit=UNIT_SUMMARY) end do ! print summary of lattices if (n_lattices > 0) then - call header("LATTICE SUMMARY") + call header("LATTICE SUMMARY", unit=UNIT_SUMMARY) do i = 1, n_lattices l => lattices(i) - call print_lattice(l) + call print_lattice(l, unit=UNIT_SUMMARY) end do end if - ! print summary of surfaces - call header("SURFACE SUMMARY") - do i = 1, n_surfaces - s => surfaces(i) - call print_surface(s) - end do - end subroutine print_geometry !=============================================================================== @@ -677,65 +724,63 @@ contains type(TallyObject), pointer :: t => null() ! Display problem summary - call header("PROBLEM SUMMARY") + call header("PROBLEM SUMMARY", unit=UNIT_SUMMARY) if (problem_type == PROB_CRITICALITY) then - write(ou,100) 'Problem type:', 'Criticality' - write(ou,101) 'Number of Cycles:', n_cycles - write(ou,101) 'Number of Inactive Cycles:', n_inactive + write(UNIT_SUMMARY,100) 'Problem type:', 'Criticality' + write(UNIT_SUMMARY,101) 'Number of Cycles:', n_cycles + write(UNIT_SUMMARY,101) 'Number of Inactive Cycles:', n_inactive elseif (problem_type == PROB_SOURCE) then - write(ou,100) 'Problem type:', 'External Source' + write(UNIT_SUMMARY,100) 'Problem type:', 'External Source' end if - write(ou,101) 'Number of Particles:', n_particles + write(UNIT_SUMMARY,101) 'Number of Particles:', n_particles ! Display geometry summary - call header("GEOMETRY SUMMARY") - write(ou,101) 'Number of Cells:', n_cells - write(ou,101) 'Number of Surfaces:', n_surfaces - write(ou,101) 'Number of Materials:', n_materials + call header("GEOMETRY SUMMARY", unit=UNIT_SUMMARY) + write(UNIT_SUMMARY,101) 'Number of Cells:', n_cells + write(UNIT_SUMMARY,101) 'Number of Surfaces:', n_surfaces + write(UNIT_SUMMARY,101) 'Number of Materials:', n_materials ! print summary of all geometry call print_geometry() ! print summary of materials - call header("MATERIAL SUMMARY") + call header("MATERIAL SUMMARY", unit=UNIT_SUMMARY) do i = 1, n_materials m => materials(i) - call print_material(m) + call print_material(m, unit=UNIT_SUMMARY) end do ! print summary of tallies if (n_tallies > 0) then - call header("TALLY SUMMARY") + call header("TALLY SUMMARY", unit=UNIT_SUMMARY) do i = 1, n_tallies t=> tallies(i) - call print_tally(t) + call print_tally(t, unit=UNIT_SUMMARY) end do end if ! print summary of unionized energy grid - call header("UNIONIZED ENERGY GRID") - write(ou,*) "Points on energy grid: " // trim(int_to_str(n_grid)) - write(ou,*) "Extra storage required: " // trim(int_to_str(& + call header("UNIONIZED ENERGY GRID", unit=UNIT_SUMMARY) + write(UNIT_SUMMARY,*) "Points on energy grid: " // trim(int_to_str(n_grid)) + write(UNIT_SUMMARY,*) "Extra storage required: " // trim(int_to_str(& n_grid*n_nuclides_total*4)) // " bytes" ! print summary of variance reduction - call header("VARIANCE REDUCTION") + call header("VARIANCE REDUCTION", unit=UNIT_SUMMARY) if (survival_biasing) then - write(ou,100) "Survival Biasing:", "on" + write(UNIT_SUMMARY,100) "Survival Biasing:", "on" else - write(ou,100) "Survival Biasing:", "off" + write(UNIT_SUMMARY,100) "Survival Biasing:", "off" end if string = real_to_str(weight_cutoff) - write(ou,100) "Weight Cutoff:", trim(string) + write(UNIT_SUMMARY,100) "Weight Cutoff:", trim(string) string = real_to_str(weight_survive) - write(ou,100) "Survival weight:", trim(string) - write(ou,*) + write(UNIT_SUMMARY,100) "Survival weight:", trim(string) ! Format descriptor for columns 100 format (1X,A,T35,A) 101 format (1X,A,T35,I11) - end subroutine print_summary !=============================================================================== @@ -811,4 +856,30 @@ contains end subroutine print_runtime +!=============================================================================== +! CREATE_SUMMARY_FILE opens the summary.out file for logging information about +! the simulation +!=============================================================================== + + subroutine create_summary_file() + + integer :: io_error + logical :: file_exists ! does log file already exist? + character(MAX_FILE_LEN) :: path ! path of summary file + + ! Create filename for log file + path = "summary.out" + + ! Check if log file already exists + inquire(FILE=path, EXIST=file_exists) + if (file_exists) then + ! Possibly copy old log file + end if + + ! Open log file for writing + open(UNIT=UNIT_SUMMARY, FILE=path, STATUS='replace', & + ACTION='write', IOSTAT=io_error) + + end subroutine create_summary_file + end module output From 409b30cf44329db852bc4c04a8f13014d81e6455 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 12 Dec 2011 09:36:35 -0600 Subject: [PATCH 02/33] Fixed typo in output module. --- src/output.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/output.F90 b/src/output.F90 index 371f239474..68ddee6bc5 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -58,7 +58,7 @@ contains ! Write information to summary file call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1) write(UNIT=UNIT_SUMMARY, FMT=*) & - "Copyright: 2010 Massachusetts Institute of Technology" + "Copyright: 2011 Massachusetts Institute of Technology" write(UNIT=UNIT_SUMMARY, FMT='(1X,A,5X,2(I1,"."),I1)') & "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",3X,A,1X,A)') & From ed4c0ef7df066133584510484fad932dbd8cad0d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 13 Dec 2011 12:18:15 -0600 Subject: [PATCH 03/33] Fixed typo in random_lcg module. --- src/random_lcg.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index 8b7225c7ff..674be52122 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -23,7 +23,7 @@ module random_lcg contains !=============================================================================== -! PRN generates a pseudo-random number using a linean congruential generator +! PRN generates a pseudo-random number using a linear congruential generator !=============================================================================== function prn() result(pseudo_rn) From 412bfe0ead34503d5007e3fd4284bb7728451617 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 13 Dec 2011 15:14:50 -0600 Subject: [PATCH 04/33] Added summary output for cross section tables and fixed small error in initialization for reactions. --- src/cross_section.F90 | 19 ++++++++--- src/output.F90 | 79 ++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 92 insertions(+), 6 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 86d15ef1d1..1ab4e41d51 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -12,7 +12,7 @@ module cross_section use fission, only: nu_total use global use material_header, only: Material - use output, only: write_message + use output, only: write_message, print_nuclide, header use string, only: split_string, str_to_int, str_to_real, & lower_case, int_to_str @@ -137,6 +137,9 @@ contains ! ========================================================================== ! READ ALL ACE CROSS SECTION TABLES + ! display header in summary.out + if (master) call header("CROSS SECTION TABLES", unit=UNIT_SUMMARY) + call dict_create(already_read) ! Loop over all files @@ -154,6 +157,10 @@ contains call read_ace_table(index_nuclides, index_list) + ! Print out information on table to summary.out file + nuc => nuclides(index_nuclides) + if (master) call print_nuclide(nuc, unit=UNIT_SUMMARY) + call dict_add_key(already_read, name, 0) call dict_add_key(already_read, alias, 0) end if @@ -634,6 +641,8 @@ contains rxn % Q_value = ZERO rxn % TY = 1 rxn % IE = 1 + rxn % has_angle_dist = .false. + rxn % has_energy_dist = .false. allocate(rxn % sigma(nuc % n_grid)) rxn % sigma = nuc % elastic @@ -650,6 +659,10 @@ contains do i = 1, NMT rxn => nuc % reactions(i+1) + ! set defaults + rxn % has_angle_dist = .false. + rxn % has_energy_dist = .false. + ! read MT number, Q-value, and neutrons produced rxn % MT = int(XSS(LMT + i - 1)) rxn % Q_value = XSS(JXS4 + i - 1) @@ -701,10 +714,6 @@ contains nuc % index_fission(i_fission) = i + 1 nuc % n_fission = nuc % n_fission + 1 end if - - ! set defaults - rxn % has_angle_dist = .false. - rxn % has_energy_dist = .false. end do end subroutine read_reactions diff --git a/src/output.F90 b/src/output.F90 index 68ddee6bc5..e5e0c5facb 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -3,7 +3,7 @@ module output use ISO_FORTRAN_ENV use constants - use cross_section_header, only: Reaction + use cross_section_header, only: Nuclide, Reaction use datatypes, only: dict_get_key use endf, only: reaction_name use geometry_header, only: Cell, Universe, Surface @@ -711,6 +711,83 @@ contains end subroutine print_geometry +!=============================================================================== +! PRINT_NUCLIDE displays information about a continuous-energy neutron +! cross_section table and its reactions and secondary angle/energy distributions +!=============================================================================== + + subroutine print_nuclide(nuc, unit) + + type(Nuclide), pointer :: nuc + integer, optional :: unit + + integer :: i + integer :: unit_ + integer :: size_total + integer :: size_xs + integer :: size_angle + integer :: size_energy + type(Reaction), pointer :: rxn => null() + + ! set default unit for writing information + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + ! Determine size of cross-sections + size_xs = (5 + nuc % n_reaction) * nuc % n_grid * 8 + size_total = size_xs + + ! Basic nuclide information + write(unit_,*) 'Nuclide ' // trim(nuc % name) + write(unit_,*) ' zaid = ' // trim(int_to_str(nuc % zaid)) + write(unit_,*) ' awr = ' // trim(real_to_str(nuc % awr)) + write(unit_,*) ' kT = ' // trim(real_to_str(nuc % kT)) + write(unit_,*) ' # of grid points = ' // trim(int_to_str(nuc % n_grid)) + write(unit_,*) ' Fissionable = ', nuc % fissionable + write(unit_,*) ' # of fission reactions = ' // trim(int_to_str(nuc % n_fission)) + write(unit_,*) ' # of reactions = ' // trim(int_to_str(nuc % n_reaction)) + write(unit_,*) ' Size of cross sections = ' // trim(int_to_str(& + size_xs)) // ' bytes' + + write(unit_,*) ' Reaction Q-value Mult IE size(angle) size(energy)' + do i = 1, nuc % n_reaction + ! Information on each reaction + rxn => nuc % reactions(i) + + ! Determine size of angle distribution + if (rxn % has_angle_dist) then + size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8 + else + size_angle = 0 + end if + + ! Determine size of energy distribution + if (rxn % has_energy_dist) then + size_energy = size(rxn % edist % data) * 8 + else + size_energy = 0 + end if + + write(unit_,'(3X,A11,1X,F8.3,2X,I3,3X,I6,1X,I11,1X,I11)') & + reaction_name(rxn % MT), rxn % Q_value, rxn % TY, rxn % IE, & + size_angle, size_energy + + ! Accumulate data size + size_total = size_total + size_angle + size_energy + end do + + ! Write total memory used + write(unit_,*) ' Total memory used = ' // trim(int_to_str(size_total)) & + // ' bytes' + + ! Blank line at end of nuclide + write(unit_,*) + + end subroutine print_nuclide + !=============================================================================== ! PRINT_SUMMARY displays summary information about the problem about to be run ! after reading all input files From 4ddcf54387e80a8b9795a3624ad9e3eae6e1c484 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 13 Dec 2011 15:25:04 -0600 Subject: [PATCH 05/33] Display number of MPI processes on stdout and summary.out. --- src/output.F90 | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/output.F90 b/src/output.F90 index e5e0c5facb..765feb2aa7 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -52,18 +52,26 @@ contains ! Write the date and time call get_today(today_date, today_time) - write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A,1X,A/)') & + write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A,1X,A)') & trim(today_date), trim(today_time) ! Write information to summary file call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1) write(UNIT=UNIT_SUMMARY, FMT=*) & - "Copyright: 2011 Massachusetts Institute of Technology" - write(UNIT=UNIT_SUMMARY, FMT='(1X,A,5X,2(I1,"."),I1)') & + "Copyright: 2011 Massachusetts Institute of Technology" + write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') & "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE - write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",3X,A,1X,A)') & + write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",5X,A,1X,A)') & trim(today_date), trim(today_time) + ! Write information on number of processors +#ifdef MPI + write(UNIT=OUTPUT_UNIT, FMT='(1X,A)') ' MPI Processes: ' // & + trim(int_to_str(n_procs)) + write(UNIT=UNIT_SUMMARY, FMT='(1X,"MPI Processes:",1X,A)') & + trim(int_to_str(n_procs)) +#endif + end subroutine title !=============================================================================== From 0617d50817fb71e2304585aef91053aca3fad842 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 13 Dec 2011 18:15:17 -0600 Subject: [PATCH 06/33] Fixed display of header for tally output. --- src/tally.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tally.F90 b/src/tally.F90 index 547f212690..29f144ac3e 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1043,7 +1043,7 @@ contains t => tallies(i) ! Write header block - call header("TALLY " // trim(int_to_str(t % id)), 3, UNIT_TALLY) + call header("TALLY " // trim(int_to_str(t % id)), unit=UNIT_TALLY, level=3) ! Handle surface current tallies separately if (t % surface_current) then From de24c7e8edb11e7cc8979331f3122b0a6abb17bf Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Dec 2011 00:00:33 -0600 Subject: [PATCH 07/33] Initial routine for calculating Shannon entropy. --- src/DEPENDENCIES | 4 +++ src/OBJECTS | 1 + src/global.F90 | 7 ++++ src/intercycle.F90 | 80 ++++++++++++++++++++++++++++++++++++++++++++++ src/main.F90 | 4 +++ 5 files changed, 96 insertions(+) create mode 100644 src/intercycle.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index f6585ba4dc..7b80cc9a0b 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -95,6 +95,9 @@ input_xml.o: xml-fortran/templates/materials_t.o input_xml.o: xml-fortran/templates/settings_t.o input_xml.o: xml-fortran/templates/tallies_t.o +intercycle.o: global.o +intercycle.o: error.o + interpolation.o: constants.o interpolation.o: endf_header.o interpolation.o: error.o @@ -105,6 +108,7 @@ interpolation.o: string.o main.o: constants.o main.o: global.o main.o: initialize.o +main.o: intercycle.o main.o: mpi_routines.o main.o: output.o main.o: particle_header.o diff --git a/src/OBJECTS b/src/OBJECTS index e8df9f0f7e..09f29ba6f8 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -15,6 +15,7 @@ geometry.o \ geometry_header.o \ global.o \ initialize.o \ +intercycle.o \ interpolation.o \ input_xml.o \ main.o \ diff --git a/src/global.F90 b/src/global.F90 index 814313c71e..4cfb6a6758 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -111,6 +111,13 @@ module global real(8) :: keff = ONE real(8) :: keff_std + ! Shannon entropy + logical :: entropy_on = .false. + real(8) :: entropy ! value of shannon entropy + real(8) :: entropy_lower_left(3) ! lower-left corner for entropy box + real(8) :: entropy_upper_right(3) ! upper-right corner for entropy box + real(8), allocatable :: entropy_p(:,:,:) + ! ============================================================================ ! PARALLEL PROCESSING VARIABLES diff --git a/src/intercycle.F90 b/src/intercycle.F90 new file mode 100644 index 0000000000..9b8cf25a66 --- /dev/null +++ b/src/intercycle.F90 @@ -0,0 +1,80 @@ +module intercycle + + use global + use error, only: warning + +#ifdef MPI + use mpi +#endif + +contains + +!=============================================================================== +! SHANNON_ENTROPY calculates the Shannon entropy of the fission source +! distribution to assess source convergence +!=============================================================================== + + subroutine shannon_entropy() + + integer :: i_bank ! index for bank sites + integer :: i ! x-index for entropy mesh + integer :: j ! y-index for entropy mesh + integer :: k ! z-index for entropy mesh + integer, save :: n_box ! total # of boxes on mesh + integer, save :: n ! # of boxes in each dimension + real(8), save :: width(3) ! width of box in each dimension + + ! On the first pass through this subroutine, we need to determine how big + ! the entropy mesh should be in each direction and then allocate a + ! three-dimensional array to store the fraction of source sites in each mesh + ! box + + if (.not. allocated(entropy_p)) then + ! determine number of boxes in each direction + n = ceiling((n_particles/20)**(1.0/3.0)) + n_box = n*n*n + + ! determine width + width = (entropy_upper_right - entropy_lower_left)/n + + ! allocate p + allocate(entropy_p(n,n,n)) + end if + + ! initialize p + entropy_p = ZERO + + ! loop over fission sites and count how many are in each mesh box + FISSION_SITES: do i_bank = 1, n_bank + ! determine indices for entropy mesh box + i = (fission_bank(i_bank) % xyz(1) - entropy_lower_left(1))/n + j = (fission_bank(i_bank) % xyz(2) - entropy_lower_left(2))/n + k = (fission_bank(i_bank) % xyz(3) - entropy_lower_left(3))/n + + ! if outside mesh, skip particle + if (i < 1 .or. i > n .or. j < 1 .or. & + j > n .or. k < 1 .or. k > n) then + message = "Fission source site outside of entropy box." + call warning() + cycle + end if + + ! add to appropriate mesh box + entropy_p(i,j,k) = entropy_p(i,j,k) + 1 + end do FISSION_SITES + + ! normalize to number of fission sites + entropy_p = entropy_p/n + + ! collect values from all processors +#ifdef MPI + call MPI_REDUCE(MPI_IN_PLACE, entropy_p, n_box, MPI_REAL8, MPI_SUM, & + 0, MPI_COMM_WORLD, mpi_err) +#endif + + ! sum values to obtain shannon entropy + if (master) entropy = sum(entropy_p * log(entropy_p)/log(2.0)) + + end subroutine shannon_entropy + +end module intercycle diff --git a/src/main.F90 b/src/main.F90 index 35479e9ea1..d7b11eb54e 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -3,6 +3,7 @@ program main use constants use global use initialize, only: initialize_run + use intercycle, only: shannon_entropy use mpi_routines, only: synchronize_bank use output, only: write_message, header, print_runtime use particle_header, only: Particle @@ -117,6 +118,9 @@ contains call timer_stop(time_ic_tallies) end if + ! Calculate shannon entropy + if (entropy_on) call shannon_entropy() + ! Distribute fission bank across processors evenly call synchronize_bank(i_cycle) From 69f8f9bd4632ec9c675074a1f5c780bca1dc8e69 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Dec 2011 08:09:36 -0600 Subject: [PATCH 08/33] Added user input option entropy_box and fixed bug in shannon_entropy subroutine. --- src/input_xml.F90 | 24 ++++++++++++++++++++++++ src/intercycle.F90 | 21 ++++++++++++++------- src/xml-fortran/templates/settings_t.xml | 1 + 3 files changed, 39 insertions(+), 7 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 078b00bc76..e05edd524a 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -142,6 +142,30 @@ contains trace_particle = trace_(2) end if + ! Entropy box + if (associated(entropy_box_)) then + ! Check to make sure enough values were supplied + if (size(entropy_box_) /= 6) then + message = "Need to supply lower-left and upper-right coordinates " & + // "for Shannon entropy box." + call fatal_error() + end if + + ! Copy values + entropy_lower_left = entropy_box_(1:3) + entropy_upper_right = entropy_box_(4:6) + + ! Check on values provided + if (.not. all(entropy_upper_right > entropy_lower_left)) then + message = "Upper-right coordinate must be greater than lower-left " & + // "coordinate for Shannon entropy box." + call fatal_error() + end if + + ! Turn on Shannon entropy calculation + entropy_on = .true. + end if + end subroutine read_settings_xml !=============================================================================== diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 9b8cf25a66..74b728768e 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -16,13 +16,14 @@ contains subroutine shannon_entropy() - integer :: i_bank ! index for bank sites + integer :: m ! index for bank sites integer :: i ! x-index for entropy mesh integer :: j ! y-index for entropy mesh integer :: k ! z-index for entropy mesh integer, save :: n_box ! total # of boxes on mesh integer, save :: n ! # of boxes in each dimension real(8), save :: width(3) ! width of box in each dimension + logical :: outside_box ! were there sites outside entropy box? ! On the first pass through this subroutine, we need to determine how big ! the entropy mesh should be in each direction and then allocate a @@ -43,19 +44,19 @@ contains ! initialize p entropy_p = ZERO + outside_box = .false. ! loop over fission sites and count how many are in each mesh box - FISSION_SITES: do i_bank = 1, n_bank + FISSION_SITES: do m = 1, n_bank ! determine indices for entropy mesh box - i = (fission_bank(i_bank) % xyz(1) - entropy_lower_left(1))/n - j = (fission_bank(i_bank) % xyz(2) - entropy_lower_left(2))/n - k = (fission_bank(i_bank) % xyz(3) - entropy_lower_left(3))/n + i = (fission_bank(m) % xyz(1) - entropy_lower_left(1))/width(1) + 1 + j = (fission_bank(m) % xyz(2) - entropy_lower_left(2))/width(2) + 1 + k = (fission_bank(m) % xyz(3) - entropy_lower_left(3))/width(3) + 1 ! if outside mesh, skip particle if (i < 1 .or. i > n .or. j < 1 .or. & j > n .or. k < 1 .or. k > n) then - message = "Fission source site outside of entropy box." - call warning() + outside_box = .true. cycle end if @@ -63,6 +64,12 @@ contains entropy_p(i,j,k) = entropy_p(i,j,k) + 1 end do FISSION_SITES + ! display warning message if there were sites outside entropy box + if (outside_box) then + message = "Fission source site(s) outside of entropy box." + call warning() + end if + ! normalize to number of fission sites entropy_p = entropy_p/n diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml index 71ae91bd27..c69439612e 100644 --- a/src/xml-fortran/templates/settings_t.xml +++ b/src/xml-fortran/templates/settings_t.xml @@ -26,5 +26,6 @@ + From 1da3280eba23c73789858885393ffdfbef52d9a9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Dec 2011 10:50:24 -0600 Subject: [PATCH 09/33] Finished Shannon entropy implementation. Closes github #17. --- src/DEPENDENCIES | 1 + src/intercycle.F90 | 131 ++++++++++++++++++++++++++++++++++++++++----- src/main.F90 | 12 +++-- src/tally.F90 | 79 +-------------------------- 4 files changed, 129 insertions(+), 94 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 7b80cc9a0b..8190c78b61 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -97,6 +97,7 @@ input_xml.o: xml-fortran/templates/tallies_t.o intercycle.o: global.o intercycle.o: error.o +intercycle.o: output.o interpolation.o: constants.o interpolation.o: endf_header.o diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 74b728768e..95dcd7003e 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -1,7 +1,10 @@ module intercycle + use ISO_FORTRAN_ENV + use global - use error, only: warning + use error, only: warning + use output, only: write_message #ifdef MPI use mpi @@ -16,10 +19,11 @@ contains subroutine shannon_entropy() - integer :: m ! index for bank sites integer :: i ! x-index for entropy mesh integer :: j ! y-index for entropy mesh integer :: k ! z-index for entropy mesh + integer :: m ! index for bank sites + integer(8) :: total_bank ! total # of fission bank sites integer, save :: n_box ! total # of boxes on mesh integer, save :: n ! # of boxes in each dimension real(8), save :: width(3) ! width of box in each dimension @@ -47,11 +51,11 @@ contains outside_box = .false. ! loop over fission sites and count how many are in each mesh box - FISSION_SITES: do m = 1, n_bank + FISSION_SITES: do m = 1, int(n_bank,4) ! determine indices for entropy mesh box - i = (fission_bank(m) % xyz(1) - entropy_lower_left(1))/width(1) + 1 - j = (fission_bank(m) % xyz(2) - entropy_lower_left(2))/width(2) + 1 - k = (fission_bank(m) % xyz(3) - entropy_lower_left(3))/width(3) + 1 + i = int((fission_bank(m) % xyz(1) - entropy_lower_left(1))/width(1)) + 1 + j = int((fission_bank(m) % xyz(2) - entropy_lower_left(2))/width(2)) + 1 + k = int((fission_bank(m) % xyz(3) - entropy_lower_left(3))/width(3)) + 1 ! if outside mesh, skip particle if (i < 1 .or. i > n .or. j < 1 .or. & @@ -70,18 +74,119 @@ contains call warning() end if - ! normalize to number of fission sites - entropy_p = entropy_p/n - - ! collect values from all processors #ifdef MPI - call MPI_REDUCE(MPI_IN_PLACE, entropy_p, n_box, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + ! collect values from all processors + if (master) then + call MPI_REDUCE(MPI_IN_PLACE, entropy_p, n_box, MPI_REAL8, MPI_SUM, & + 0, MPI_COMM_WORLD, mpi_err) + else + call MPI_REDUCE(entropy_p, entropy_p, n_box, MPI_REAL8, MPI_SUM, & + 0, MPI_COMM_WORLD, mpi_err) + end if + + ! determine total number of bank sites + call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) +#else + total_bank = n_bank #endif ! sum values to obtain shannon entropy - if (master) entropy = sum(entropy_p * log(entropy_p)/log(2.0)) + if (master) then + entropy_p = entropy_p / total_bank + entropy = -sum(entropy_p * log(entropy_p)/log(2.0), entropy_p > ZERO) + end if end subroutine shannon_entropy +!=============================================================================== +! CALCULATE_KEFF calculates the single cycle estimate of keff as well as the +! mean and standard deviation of the mean for active cycles and displays them +!=============================================================================== + + subroutine calculate_keff(i_cycle) + + integer, intent(in) :: i_cycle ! index of current cycle + + integer(8) :: total_bank ! total number of source sites + integer :: n ! active cycle number + real(8) :: k_cycle ! single cycle estimate of keff + real(8), save :: k_sum ! accumulated keff + real(8), save :: k_sum_sq ! accumulated keff**2 + + message = "Calculate cycle keff..." + call write_message(8) + + ! initialize sum and square of sum at beginning of run + if (i_cycle == 1) then + k_sum = ZERO + k_sum_sq = ZERO + end if + +#ifdef MPI + ! Collect number bank sites onto master process + call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & + & MPI_COMM_WORLD, mpi_err) +#else + total_bank = n_bank +#endif + + ! Collect statistics and print output + if (master) then + ! Since the creation of bank sites was originally weighted by the last + ! cycle keff, we need to multiply by that keff to get the current cycle's + ! value + + k_cycle = real(total_bank)/real(n_particles)*keff + + if (i_cycle > n_inactive) then + ! Active cycle number + n = i_cycle - n_inactive + + ! Accumulate cycle estimate of k + k_sum = k_sum + k_cycle + k_sum_sq = k_sum_sq + k_cycle*k_cycle + + ! Determine mean and standard deviation of mean + keff = k_sum/n + keff_std = sqrt((k_sum_sq/n - keff*keff)/n) + + ! Display output for this cycle + if (i_cycle > n_inactive+1) then + if (entropy_on) then + write(UNIT=OUTPUT_UNIT, FMT=103) i_cycle, k_cycle, entropy, & + keff, keff_std + else + write(UNIT=OUTPUT_UNIT, FMT=101) i_cycle, k_cycle, keff, keff_std + end if + else + if (entropy_on) then + write(UNIT=OUTPUT_UNIT, FMT=102) i_cycle, k_cycle, entropy + else + write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle + end if + end if + else + ! Display output for inactive cycle + if (entropy_on) then + write(UNIT=OUTPUT_UNIT, FMT=102) i_cycle, k_cycle, entropy + else + write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle + end if + keff = k_cycle + end if + end if + +#ifdef MPI + ! Broadcast new keff value to all processors + call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err) +#endif + +100 format (2X,I5,2X,F8.5) +101 format (2X,I5,2X,F8.5,9X,F8.5," +/-",F8.5) +102 format (2X,I5,2X,F8.5,3X,F8.5) +103 format (2X,I5,2X,F8.5,3X,F8.5,3X,F8.5," +/-",F8.5) + + end subroutine calculate_keff + end module intercycle diff --git a/src/main.F90 b/src/main.F90 index d7b11eb54e..c9534c09e1 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -3,7 +3,7 @@ program main use constants use global use initialize, only: initialize_run - use intercycle, only: shannon_entropy + use intercycle, only: shannon_entropy, calculate_keff use mpi_routines, only: synchronize_bank use output, only: write_message, header, print_runtime use particle_header, only: Particle @@ -13,7 +13,7 @@ program main use source, only: get_source_particle use string, only: int_to_str use tally, only: synchronize_tallies, write_tallies, & - tally_statistics, calculate_keff + tally_statistics use timing, only: timer_start, timer_stop #ifdef MPI @@ -64,6 +64,12 @@ contains tallies_on = .false. call timer_start(time_inactive) + ! Display column titles + message = " Cycle k(cycle) Entropy Average k" + call write_message(1) + message = " ===== ======== ======= ===================" + call write_message(1) + ! ========================================================================== ! LOOP OVER CYCLES CYCLE_LOOP: do i_cycle = 1, n_cycles @@ -73,7 +79,7 @@ contains message = "Simulating cycle " // trim(int_to_str(i_cycle)) // "..." call write_message(8) - + ! Set all tallies to zero n_bank = 0 diff --git a/src/tally.F90 b/src/tally.F90 index 29f144ac3e..8b039c8af6 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1,13 +1,11 @@ module tally - use ISO_FORTRAN_ENV - use constants use error, only: fatal_error use global use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices use mesh_header, only: StructuredMesh - use output, only: write_message, header + use output, only: header use search, only: binary_search use string, only: int_to_str, real_to_str use tally_header, only: TallyScore, TallyMapItem, TallyMapElement @@ -23,81 +21,6 @@ module tally contains -!=============================================================================== -! CALCULATE_KEFF calculates the single cycle estimate of keff as well as the -! mean and standard deviation of the mean for active cycles and displays them -!=============================================================================== - - subroutine calculate_keff(i_cycle) - - integer, intent(in) :: i_cycle ! index of current cycle - - integer(8) :: total_bank ! total number of source sites - integer :: n ! active cycle number - real(8) :: k_cycle ! single cycle estimate of keff - real(8), save :: k_sum ! accumulated keff - real(8), save :: k_sum_sq ! accumulated keff**2 - - message = "Calculate cycle keff..." - call write_message(8) - - ! initialize sum and square of sum at beginning of run - if (i_cycle == 1) then - k_sum = ZERO - k_sum_sq = ZERO - end if - -#ifdef MPI - ! Collect number bank sites onto master process - call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & - & MPI_COMM_WORLD, mpi_err) -#else - total_bank = n_bank -#endif - - ! Collect statistics and print output - if (master) then - ! Since the creation of bank sites was originally weighted by the last - ! cycle keff, we need to multiply by that keff to get the current cycle's - ! value - - k_cycle = real(total_bank)/real(n_particles)*keff - - if (i_cycle > n_inactive) then - ! Active cycle number - n = i_cycle - n_inactive - - ! Accumulate cycle estimate of k - k_sum = k_sum + k_cycle - k_sum_sq = k_sum_sq + k_cycle*k_cycle - - ! Determine mean and standard deviation of mean - keff = k_sum/n - keff_std = sqrt((k_sum_sq/n - keff*keff)/n) - - ! Display output for this cycle - if (i_cycle > n_inactive+1) then - write(UNIT=OUTPUT_UNIT, FMT=101) i_cycle, k_cycle, keff, keff_std - else - write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle - end if - else - ! Display output for inactive cycle - write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle - keff = k_cycle - end if - end if - -#ifdef MPI - ! Broadcast new keff value to all processors - call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err) -#endif - -100 format (2X,I5,2X,F8.5) -101 format (2X,I5,2X,F8.5,9X,F8.5,1X,F8.5) - - end subroutine calculate_keff - !=============================================================================== ! CREATE_TALLY_MAP creates a map that allows a quick determination of which ! tallies and bins need to be scored to when a particle makes a collision. This From 1a7943c4c81bfcb8d6276cadcef637bd9da9d171 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Dec 2011 11:14:55 -0600 Subject: [PATCH 10/33] Fixed column titles for cases without Shannon entropy. --- src/intercycle.F90 | 2 +- src/main.F90 | 15 +++++++++++---- src/mpi_routines.F90 | 3 ++- 3 files changed, 14 insertions(+), 6 deletions(-) diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 95dcd7003e..3ec77bf6e0 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -183,7 +183,7 @@ contains #endif 100 format (2X,I5,2X,F8.5) -101 format (2X,I5,2X,F8.5,9X,F8.5," +/-",F8.5) +101 format (2X,I5,2X,F8.5,5X,F8.5," +/-",F8.5) 102 format (2X,I5,2X,F8.5,3X,F8.5) 103 format (2X,I5,2X,F8.5,3X,F8.5,3X,F8.5," +/-",F8.5) diff --git a/src/main.F90 b/src/main.F90 index c9534c09e1..ca0094b140 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -65,10 +65,17 @@ contains call timer_start(time_inactive) ! Display column titles - message = " Cycle k(cycle) Entropy Average k" - call write_message(1) - message = " ===== ======== ======= ===================" - call write_message(1) + if (entropy_on) then + message = " Cycle k(cycle) Entropy Average k" + call write_message(1) + message = " ===== ======== ======= ===================" + call write_message(1) + else + message = " Cycle k(cycle) Average k" + call write_message(1) + message = " ===== ======== ===================" + call write_message(1) + end if ! ========================================================================== ! LOOP OVER CYCLES diff --git a/src/mpi_routines.F90 b/src/mpi_routines.F90 index a523694772..3729fc98f4 100644 --- a/src/mpi_routines.F90 +++ b/src/mpi_routines.F90 @@ -386,7 +386,8 @@ contains tally_temp = t % scores(:,:) % val_history if (master) then - ! Description of MPI_IN_PLANE + ! The MPI_IN_PLACE specifier allows the master to copy values into a + ! receive buffer without having a temporary variable call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, MPI_SUM, & 0, MPI_COMM_WORLD, mpi_err) From ecca5388a1c3bafc03f1f3308650b0257b540978 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 15 Dec 2011 13:29:00 -0600 Subject: [PATCH 11/33] Fixed redundant calculation of lattice boundary distance. Reduced FP_PRECISION. --- src/constants.F90 | 9 +++++- src/geometry.F90 | 77 ++++++++++++++--------------------------------- src/physics.F90 | 9 +++--- src/plot.F90 | 6 ++-- 4 files changed, 37 insertions(+), 64 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 352d500c4b..4fdaf1725a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -22,7 +22,7 @@ module constants real(8), parameter :: TINY_BIT = 1e-8_8 ! User for precision in geometry - real(8), parameter :: FP_PRECISION = 1e-7_8 + real(8), parameter :: FP_PRECISION = 1e-5_8 ! Maximum number of collisions/crossings integer, parameter :: MAX_EVENTS = 10000 @@ -79,6 +79,13 @@ module constants & LATTICE_RECT = 1, & & LATTICE_HEX = 2 + ! Lattice boundary crossings + integer, parameter :: & + & LATTICE_LEFT = 1, & + & LATTICE_RIGHT = 2, & + & LATTICE_BOTTOM = 3, & + & LATTICE_TOP = 4 + ! Surface types integer, parameter :: & & SURF_PX = 1, & ! Plane parallel to x-plane diff --git a/src/geometry.F90 b/src/geometry.F90 index 569755d2e1..3ac71e768c 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -424,22 +424,13 @@ contains ! CROSS_LATTICE moves a particle into a new lattice element !=============================================================================== - subroutine cross_lattice(p) + subroutine cross_lattice(p, lattice_crossed) type(Particle), pointer :: p + integer, intent(in) :: lattice_crossed integer :: i_x ! x index in lattice integer :: i_y ! y index in lattice - real(8) :: d_left ! distance to left side - real(8) :: d_right ! distance to right side - real(8) :: d_bottom ! distance to bottom side - real(8) :: d_top ! distance to top side - real(8) :: dist ! shortest distance - real(8) :: x ! x coordinate in local lattice element - real(8) :: y ! y coordinate in local lattice element - real(8) :: z ! z coordinate in local lattice element - real(8) :: u ! cosine of angle with x axis - real(8) :: v ! cosine of angle with y axis real(8) :: x0 ! half the width of lattice element real(8) :: y0 ! half the height of lattice element logical :: found ! particle found in cell? @@ -454,68 +445,36 @@ contains call write_message() end if - u = p % coord % uvw(1) - v = p % coord % uvw(2) - if (lat % type == LATTICE_RECT) then - x = p % coord % xyz(1) - y = p % coord % xyz(2) - z = p % coord % xyz(3) x0 = lat % width_x * 0.5_8 y0 = lat % width_y * 0.5_8 - - dist = INFINITY - ! left and right sides - if (u == ZERO) then - d_left = INFINITY - d_right = INFINITY - elseif (u > 0) then - d_left = INFINITY - d_right = (x0 - x)/u - else - d_left = -(x0 + x)/u - d_right = INFINITY - end if - - ! top and bottom sides - if (v == ZERO) then - d_bottom = INFINITY - d_top = INFINITY - elseif (v > 0) then - d_bottom = INFINITY - d_top = (y0 - y)/v - else - d_bottom = -(y0 + y)/v - d_top = INFINITY - end if - - dist = min(d_left, d_right, d_top, d_bottom) - if (dist == d_left) then + select case (lattice_crossed) + case (LATTICE_LEFT) ! Move particle to left element p % coord % lattice_x = p % coord % lattice_x - 1 p % coord % xyz(1) = x0 - elseif (dist == d_right) then + case (LATTICE_RIGHT) ! Move particle to right element p % coord % lattice_x = p % coord % lattice_x + 1 p % coord % xyz(1) = -x0 - elseif (dist == d_bottom) then + case (LATTICE_BOTTOM) ! Move particle to bottom element p % coord % lattice_y = p % coord % lattice_y - 1 p % coord % xyz(2) = y0 - elseif (dist == d_top) then + case (LATTICE_TOP) ! Move particle to top element p % coord % lattice_y = p % coord % lattice_y + 1 p % coord % xyz(2) = -y0 - end if + end select elseif (lat % type == LATTICE_HEX) then ! TODO: Add hex lattice support end if - + ! Check to make sure still in lattice i_x = p % coord % lattice_x i_y = p % coord % lattice_y @@ -555,7 +514,7 @@ contains type(Particle), pointer :: p real(8), intent(out) :: dist integer, intent(out) :: surface_crossed - logical, intent(out) :: lattice_crossed + integer, intent(out) :: lattice_crossed integer :: i ! index for surface in cell integer :: index_surf ! index in surfaces array (with sign) @@ -576,7 +535,7 @@ contains ! inialize distance to infinity (huge) dist = INFINITY - lattice_crossed = .false. + lattice_crossed = NONE nullify(final_coord) ! Get pointer to top-level coordinates @@ -866,7 +825,7 @@ contains if (d < dist) then dist = d surface_crossed = -cl % surfaces(i) - lattice_crossed = .false. + lattice_crossed = NONE final_coord => coord end if @@ -906,7 +865,11 @@ contains if (d < dist) then if (abs(d - dist)/dist >= FP_PRECISION) then dist = d - lattice_crossed = .true. + if (u > 0) then + lattice_crossed = LATTICE_RIGHT + else + lattice_crossed = LATTICE_LEFT + end if final_coord => coord end if end if @@ -923,7 +886,11 @@ contains if (d < dist) then if (abs(d - dist)/dist >= FP_PRECISION) then dist = d - lattice_crossed = .true. + if (v > 0) then + lattice_crossed = LATTICE_TOP + else + lattice_crossed = LATTICE_BOTTOM + end if final_coord => coord end if end if diff --git a/src/physics.F90 b/src/physics.F90 index 04c829ff86..040da1cd00 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -30,13 +30,13 @@ contains type(Particle), pointer :: p integer :: surface_crossed ! surface which particle is on + integer :: lattice_crossed ! lattice boundary which particle crossed integer :: last_cell ! most recent cell particle was in integer :: n_event ! number of collisions/crossings real(8) :: d_boundary ! distance to nearest boundary real(8) :: d_collision ! sampled distance to collision real(8) :: distance ! distance particle travels logical :: found_cell ! found cell which particle is in? - logical :: lattice_crossed ! is surface crossing in lattice? type(LocalCoord), pointer :: coord => null() if (p % coord % cell == NONE) then @@ -94,9 +94,9 @@ contains if (d_collision > d_boundary) then last_cell = p % coord % cell p % coord % cell = NONE - if (lattice_crossed) then + if (lattice_crossed /= NONE) then p % surface = NONE - call cross_lattice(p) + call cross_lattice(p, lattice_crossed) else p % surface = surface_crossed call cross_surface(p, last_cell) @@ -1380,8 +1380,7 @@ contains mu = mu0 + (xi - c_k)/p0 elseif (interp == LINEAR_LINEAR) then - ! Linear-linear interpolation -- not sure how you come about the - ! formula given in the MCNP manual + ! Linear-linear interpolation p1 = rxn % adist % data(lc + NP + k+1) mu1 = rxn % adist % data(lc + k+1) diff --git a/src/plot.F90 b/src/plot.F90 index a42c2ac4ed..77ae883c03 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -23,6 +23,7 @@ contains integer :: i ! loop index integer :: surface_crossed ! surface which particle is on + integer :: lattice_crossed ! is surface crossing in lattice? integer :: last_cell ! most recent cell particle was in integer :: enter_surface ! entrance surface real(8) :: xyz(3) ! starting coordinates @@ -31,7 +32,6 @@ contains real(8) :: d ! distance to boundary real(8) :: distance ! distance particle travels logical :: found_cell ! found cell which particle is in? - logical :: lattice_crossed ! is surface crossing in lattice? character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file type(Cell), pointer :: c => null() type(Universe), pointer :: univ => null() @@ -161,9 +161,9 @@ contains write(UNIT=UNIT_PLOT) p % coord0 % xyz, last_cell p % coord % cell = 0 - if (lattice_crossed) then + if (lattice_crossed /= NONE) then p % surface = NONE - call cross_lattice(p) + call cross_lattice(p, lattice_crossed) else p % surface = surface_crossed call cross_surface(p, last_cell) From 50fbdd9ea631f0e98b6eac6a3ead168eb020b769 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 15 Dec 2011 13:52:52 -0600 Subject: [PATCH 12/33] Fixed typo in sample_target_velocity. --- src/physics.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics.F90 b/src/physics.F90 index 040da1cd00..4547b9c6c7 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1023,7 +1023,7 @@ contains ! determine velocity vector of target nucleus v_target(1) = u*vt v_target(2) = v*vt - v_target(3) = W*vt + v_target(3) = w*vt end subroutine sample_target_velocity From 0431b19504527ca5e662f4e5e87d8bc6f4b9864a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 16 Dec 2011 10:28:56 -0600 Subject: [PATCH 13/33] Added out-of-bounds check on lattice elements in find_cell. --- src/geometry.F90 | 88 ++++++++++++++++++++++++++++-------------------- 1 file changed, 52 insertions(+), 36 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 3ac71e768c..29ad60626c 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -171,26 +171,42 @@ contains x = ceiling((xyz(1) - lat % x0)/lat % width_x) y = ceiling((xyz(2) - lat % y0)/lat % width_y) - ! Create new level of coordinates - p % in_lower_universe = .true. - allocate(p % coord % next) - - ! adjust local position of particle - p % coord % next % xyz(1) = p % coord % xyz(1) - & - (lat%x0 + (x-0.5_8)*lat%width_x) - p % coord % next % xyz(2) = p % coord % xyz(2) - & - (lat%y0 + (y-0.5_8)*lat%width_y) - p % coord % next % xyz(3) = p % coord % xyz(3) - p % coord % next % uvw = p % coord % uvw + ! Check if lattice coordinates are within bounds + if (x < 1 .or. x > lat % n_x .or. & + y < 1 .or. y > lat % n_y) then - ! Move particle to next level - p % coord => p % coord % next + ! This condition should only get hit in rare circumstances where + ! a neutron hits the corner of a lattice. In this case, the + ! neutron may need to be moved diagonally across the lattice. To + ! do so, we remove all lower coordinate levels and then search + ! from universe 0. + + p % coord => p % coord0 + call deallocate_coord(p % coord % next) + + else + + ! Create new level of coordinates + p % in_lower_universe = .true. + allocate(p % coord % next) - ! set particle lattice indices - p % coord % lattice = c % fill - p % coord % lattice_x = x - p % coord % lattice_y = y - p % coord % universe = lat % element(x,y) + ! adjust local position of particle + p % coord % next % xyz(1) = p % coord % xyz(1) - & + (lat%x0 + (x-0.5_8)*lat%width_x) + p % coord % next % xyz(2) = p % coord % xyz(2) - & + (lat%y0 + (y-0.5_8)*lat%width_y) + p % coord % next % xyz(3) = p % coord % xyz(3) + p % coord % next % uvw = p % coord % uvw + + ! Move particle to next level + p % coord => p % coord % next + + ! set particle lattice indices + p % coord % lattice = c % fill + p % coord % lattice_x = x + p % coord % lattice_y = y + p % coord % universe = lat % element(x,y) + end if call find_cell(p, found) if (.not. found) exit @@ -216,19 +232,19 @@ contains type(Particle), pointer :: p 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 + logical :: found ! particle found in universe? type(Surface), pointer :: surf => null() surf => surfaces(abs(p % surface)) @@ -429,11 +445,11 @@ contains type(Particle), pointer :: p integer, intent(in) :: lattice_crossed - integer :: i_x ! x index in lattice - integer :: i_y ! y index in lattice - real(8) :: x0 ! half the width of lattice element - real(8) :: y0 ! half the height of lattice element - logical :: found ! particle found in cell? + integer :: i_x ! x index in lattice + integer :: i_y ! y index in lattice + real(8) :: x0 ! half the width of lattice element + real(8) :: y0 ! half the height of lattice element + logical :: found ! particle found in cell? type(Lattice), pointer :: lat => null() lat => lattices(p % coord % lattice) From efe3feeb546c3509126386848d4753cce3b0dcd2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 16 Dec 2011 17:27:30 +0000 Subject: [PATCH 14/33] Added Makefile option for IBM Blue Gene/P. --- src/Makefile | 18 +++++++++++++----- src/xml-fortran/Makefile | 4 ++++ 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/src/Makefile b/src/Makefile index a95443ce3b..bb81589dc7 100644 --- a/src/Makefile +++ b/src/Makefile @@ -109,7 +109,7 @@ endif ifeq ($(COMPILER),ibm) F90 = xlf2003 - F90FLAGS := -WF,-DNO_F2008 + F90FLAGS := -WF,-DNO_F2008 -O2 # Debugging ifeq ($(DEBUG),yes) @@ -157,11 +157,19 @@ ifeq ($(USE_MPI),yes) endif #=============================================================================== -# Special options for ORNL Jaguar supercomputer +# Options for IBM Blue Gene/P ANL supercomputer #=============================================================================== -HOSTNAME = $(shell hostname) -ifneq (,$(findstring jaguar,$(HOSTNAME))) +ifeq ($(MACHINE),bluegene) + F90 = /bgsys/drivers/ppcfloor/comm/xl/bin/mpixlf2003 + F90FLAGS = -WF,-DNO_F2008,-DMPI -O3 +endif + +#=============================================================================== +# Options for Cray XK6 ORNL Jaguar supercomputer +#=============================================================================== + +ifeq ($(MACHINE),crayxk6) F90 = ftn F90FLAGS += -DMPI endif @@ -172,7 +180,7 @@ endif all: xml-fortran $(program) xml-fortran: - cd xml-fortran; make F90=$(F90) F90FLAGS="$(F90FLAGS)" + cd xml-fortran; make MACHINE=$(MACHINE) F90=$(F90) F90FLAGS="$(F90FLAGS)" cd xml-fortran/templates; make F90=$(F90) F90FLAGS="$(F90FLAGS)" $(program): $(objects) $(F90) $(objects) $(templates) $(xml_fort) -o $@ $(LDFLAGS) diff --git a/src/xml-fortran/Makefile b/src/xml-fortran/Makefile index 6e963f0595..c0514533a3 100644 --- a/src/xml-fortran/Makefile +++ b/src/xml-fortran/Makefile @@ -8,6 +8,10 @@ objects = $(source:.f90=.o) # Ignore unused variables +ifeq ($(MACHINE),bluegene) + override F90 = xlf2003 +endif + ifeq ($(F90),ifort) override F90FLAGS += -warn nounused endif From 56de9654adb77e4a33881a53164e0071febbf73e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 16 Dec 2011 12:54:09 -0600 Subject: [PATCH 15/33] Remove call to list_size from add_grid_points. --- src/energy_grid.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index d2351c94cf..39f5fcd31f 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -99,7 +99,7 @@ contains ! if the original list is empty, we need to allocate the first element and ! store first energy point - if (list_size(list) == 0) then + if (.not. associated(list)) then allocate(list) current => list do i = 1, n From 9e30afb7333dbc468bf48d97808a50083d10cb6d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 16 Dec 2011 16:18:53 -0600 Subject: [PATCH 16/33] Moved calculation of S(a,b) cross sections into separate subroutine. --- src/physics.F90 | 154 ++++++++++++++++++++++++++---------------------- 1 file changed, 85 insertions(+), 69 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 4547b9c6c7..0db20f28c9 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -226,13 +226,8 @@ contains integer, intent(in) :: index_sab ! index into sab_tables array integer :: IE ! index on nuclide energy grid - integer :: IE_sab ! index on S(a,b) energy grid real(8) :: f ! interp factor on nuclide energy grid - real(8) :: f_sab ! interp factor on S(a,b) energy grid - real(8) :: inelastic ! S(a,b) inelastic cross section - real(8) :: elastic ! S(a,b) elastic cross section type(Nuclide), pointer :: nuc => null() - type(SAB_Table), pointer :: sab => null() ! Set pointer to nuclide nuc => nuclides(index_nuclide) @@ -282,76 +277,97 @@ contains ! need to correct it by subtracting the non-S(a,b) elastic cross section and ! then add back in the calculated S(a,b) elastic+inelastic cross section. - if (index_sab > 0) then - micro_xs(index_nuclide) % use_sab = .true. - - ! Get pointer to S(a,b) table - sab => sab_tables(index_sab) - - ! Get index and interpolation factor for inelastic grid - if (p%E < sab % inelastic_e_in(1)) then - IE_sab = 1 - f_sab = ZERO - else - IE_sab = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E) - f_sab = (p%E - sab%inelastic_e_in(IE_sab)) / & - (sab%inelastic_e_in(IE_sab+1) - sab%inelastic_e_in(IE_sab)) - end if - - ! Calculate S(a,b) inelastic scattering cross section - inelastic = (ONE-f_sab) * sab % inelastic_sigma(IE_sab) + f_sab * & - sab % inelastic_sigma(IE_sab + 1) - - ! Check for elastic data - if (p % E < sab % threshold_elastic) then - ! Determine whether elastic scattering is given in the coherent or - ! incoherent approximation. For coherent, the cross section is - ! represented as P/E whereas for incoherent, it is simply P - - if (sab % elastic_mode == SAB_ELASTIC_EXACT) then - if (p % E < sab % elastic_e_in(1)) then - ! If energy is below that of the lowest Bragg peak, the elastic - ! cross section will be zero - elastic = ZERO - else - IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) - elastic = sab % elastic_P(IE_sab) / p % E - end if - else - ! Determine index on elastic energy grid - if (p % E < sab % elastic_e_in(1)) then - IE_sab = 1 - else - IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) - end if - - ! Get interpolation factor for elastic grid - f_sab = (p%E - sab%elastic_e_in(IE_sab))/(sab%elastic_e_in(IE_sab+1) - & - sab%elastic_e_in(IE_sab)) - - ! Calculate S(a,b) elastic scattering cross section - elastic = (ONE-f_sab) * sab % elastic_P(IE_sab) + f_sab * & - sab % elastic_P(IE_sab + 1) - end if - else - ! No elastic data - elastic = ZERO - end if - - ! Correct total and elastic cross sections - micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % total - & - micro_xs(index_nuclide) % elastic + inelastic + elastic - micro_xs(index_nuclide) % elastic = inelastic + elastic - - ! Store S(a,b) elastic cross section for sampling later - micro_xs(index_nuclide) % elastic_sab = elastic - end if + if (index_sab > 0) call calculate_sab_xs(p, index_nuclide, index_sab) ! Set last evaluated energy micro_xs(index_nuclide) % last_E = p % E end subroutine calculate_nuclide_xs +!=============================================================================== +! CALCULATE_SAB_XS determines the elastic and inelastic scattering +! cross-sections in the thermal energy range. These cross sections replace +! whatever data were taken from the normal Nuclide table. +!=============================================================================== + + subroutine calculate_sab_xs(p, index_nuclide, index_sab) + + type(Particle), pointer :: p + integer, intent(in) :: index_nuclide ! index into nuclides array + integer, intent(in) :: index_sab ! index into sab_tables array + + integer :: IE_sab ! index on S(a,b) energy grid + real(8) :: f_sab ! interp factor on S(a,b) energy grid + real(8) :: inelastic ! S(a,b) inelastic cross section + real(8) :: elastic ! S(a,b) elastic cross section + type(SAB_Table), pointer :: sab => null() + + ! Set flag that S(a,b) treatment should be used for scattering + micro_xs(index_nuclide) % use_sab = .true. + + ! Get pointer to S(a,b) table + sab => sab_tables(index_sab) + + ! Get index and interpolation factor for inelastic grid + if (p%E < sab % inelastic_e_in(1)) then + IE_sab = 1 + f_sab = ZERO + else + IE_sab = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E) + f_sab = (p%E - sab%inelastic_e_in(IE_sab)) / & + (sab%inelastic_e_in(IE_sab+1) - sab%inelastic_e_in(IE_sab)) + end if + + ! Calculate S(a,b) inelastic scattering cross section + inelastic = (ONE-f_sab) * sab % inelastic_sigma(IE_sab) + f_sab * & + sab % inelastic_sigma(IE_sab + 1) + + ! Check for elastic data + if (p % E < sab % threshold_elastic) then + ! Determine whether elastic scattering is given in the coherent or + ! incoherent approximation. For coherent, the cross section is + ! represented as P/E whereas for incoherent, it is simply P + + if (sab % elastic_mode == SAB_ELASTIC_EXACT) then + if (p % E < sab % elastic_e_in(1)) then + ! If energy is below that of the lowest Bragg peak, the elastic + ! cross section will be zero + elastic = ZERO + else + IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) + elastic = sab % elastic_P(IE_sab) / p % E + end if + else + ! Determine index on elastic energy grid + if (p % E < sab % elastic_e_in(1)) then + IE_sab = 1 + else + IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) + end if + + ! Get interpolation factor for elastic grid + f_sab = (p%E - sab%elastic_e_in(IE_sab))/(sab%elastic_e_in(IE_sab+1) - & + sab%elastic_e_in(IE_sab)) + + ! Calculate S(a,b) elastic scattering cross section + elastic = (ONE-f_sab) * sab % elastic_P(IE_sab) + f_sab * & + sab % elastic_P(IE_sab + 1) + end if + else + ! No elastic data + elastic = ZERO + end if + + ! Correct total and elastic cross sections + micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % total - & + micro_xs(index_nuclide) % elastic + inelastic + elastic + micro_xs(index_nuclide) % elastic = inelastic + elastic + + ! Store S(a,b) elastic cross section for sampling later + micro_xs(index_nuclide) % elastic_sab = elastic + + end subroutine calculate_sab_xs + !=============================================================================== ! FIND_ENERGY_INDEX determines the index on the union energy grid and the ! interpolation factor for a particle at a certain energy From d69dee4f7f0336af7441e5cc8b93ce963a6e1c5c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 16 Dec 2011 17:46:01 -0600 Subject: [PATCH 17/33] Change structure of UrrData. --- src/cross_section.F90 | 21 ++++++++++++--------- src/cross_section_header.F90 | 9 +++++++-- 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 1ab4e41d51..fde72e51b9 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -1089,7 +1089,6 @@ contains if (JXS23 /= 0) then nuc % urr_present = .true. allocate(nuc % urr_data) - allocate(nuc % urr_data % params(6)) lc = JXS23 else nuc % urr_present = .false. @@ -1097,16 +1096,20 @@ contains end if ! read parameters - nuc % urr_data % params(1) = int(XSS(lc)) ! # of incident energies - nuc % urr_data % params(2) = int(XSS(lc + 1)) ! # of probabilities - nuc % urr_data % params(3) = int(XSS(lc + 2)) ! interpolation parameter - nuc % urr_data % params(4) = int(XSS(lc + 3)) ! inelastic competition flag - nuc % urr_data % params(5) = int(XSS(lc + 4)) ! other absorption flag - nuc % urr_data % params(6) = int(XSS(lc + 5)) ! factors flag + nuc % urr_data % n_energy = int(XSS(lc)) + nuc % urr_data % n_prob = int(XSS(lc + 1)) + nuc % urr_data % interp = int(XSS(lc + 2)) + nuc % urr_data % inelastic_flag = int(XSS(lc + 3)) + nuc % urr_data % absorption_flag = int(XSS(lc + 4)) + if (int(XSS(lc + 5)) == 0) then + nuc % urr_data % multiply_smooth = .false. + else + nuc % urr_data % multiply_smooth = .true. + end if ! allocate incident energies and probability tables - N = nuc % urr_data % params(1) - M = nuc % urr_data % params(2) + N = nuc % urr_data % n_energy + M = nuc % urr_data % n_prob allocate(nuc % urr_data % energy(N)) allocate(nuc % urr_data % prob(N,6,M)) diff --git a/src/cross_section_header.F90 b/src/cross_section_header.F90 index f2a64e306f..3eaf6d4462 100644 --- a/src/cross_section_header.F90 +++ b/src/cross_section_header.F90 @@ -51,11 +51,16 @@ module cross_section_header end type Reaction !=============================================================================== -! URRDATA contains unresolved resonance data. +! URRDATA contains probability tables for the unresolved resonance range. !=============================================================================== type UrrData - integer, allocatable :: params(:) + integer :: n_energy ! # of incident neutron energies + integer :: n_prob ! # of probabilities + integer :: interp ! inteprolation (2=lin-lin, 5=log-log) + integer :: inelastic_flag ! inelastic competition flag + integer :: absorption_flag ! other absorption flag + logical :: multiply_smooth ! multiply by smooth cross section? real(8), allocatable :: energy(:) real(8), allocatable :: prob(:,:,:) end type UrrData From cf7f794381bf7adbd42c8f3e50d048bc5d7784e7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 16 Dec 2011 18:14:09 -0600 Subject: [PATCH 18/33] Added summary output for URR probablity tables. --- src/output.F90 | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/output.F90 b/src/output.F90 index 765feb2aa7..d51f83ee09 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -3,7 +3,7 @@ module output use ISO_FORTRAN_ENV use constants - use cross_section_header, only: Nuclide, Reaction + use cross_section_header, only: Nuclide, Reaction, UrrData use datatypes, only: dict_get_key use endf, only: reaction_name use geometry_header, only: Cell, Universe, Surface @@ -736,6 +736,7 @@ contains integer :: size_angle integer :: size_energy type(Reaction), pointer :: rxn => null() + type(UrrData), pointer :: urr => null() ! set default unit for writing information if (present(unit)) then @@ -787,6 +788,20 @@ contains size_total = size_total + size_angle + size_energy end do + ! Write information about URR probability tables + if (nuc % urr_present) then + urr => nuc % urr_data + write(unit_,*) ' Unresolved resonance probability table:' + write(unit_,*) ' # of energies = ' // trim(int_to_str(urr % n_energy)) + write(unit_,*) ' # of probabilities = ' // trim(int_to_str(urr % n_prob)) + write(unit_,*) ' Interpolation = ' // trim(int_to_str(urr % interp)) + write(unit_,*) ' Inelastic flag = ' // trim(int_to_str(urr % inelastic_flag)) + write(unit_,*) ' Absorption flag = ' // trim(int_to_str(urr % absorption_flag)) + write(unit_,*) ' Multiply by smooth? ', urr % multiply_smooth + write(unit_,*) ' Min energy = ', trim(real_to_str(urr % energy(1))) + write(unit_,*) ' Max energy = ', trim(real_to_str(urr % energy(urr % n_energy))) + end if + ! Write total memory used write(unit_,*) ' Total memory used = ' // trim(int_to_str(size_total)) & // ' bytes' From f31c7e5f0a795c8d3956dbabccde5684fa49fed4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 18 Dec 2011 15:47:23 -0500 Subject: [PATCH 19/33] Added most logic for unresolved resonance probability tables. Does not work yet. --- src/constants.F90 | 9 +++ src/cross_section.F90 | 19 ++++++ src/cross_section_header.F90 | 1 + src/global.F90 | 3 + src/physics.F90 | 119 ++++++++++++++++++++++++++++++++++- 5 files changed, 150 insertions(+), 1 deletion(-) diff --git a/src/constants.F90 b/src/constants.F90 index 4fdaf1725a..5188035f83 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -208,6 +208,15 @@ module constants ASCII = 1, & ! ASCII cross section file BINARY = 2 ! Binary cross section file + ! Probability table parameters + integer, parameter :: & + URR_CUM_PROB = 1, & + URR_TOTAL = 2, & + URR_ELASTIC = 3, & + URR_FISSION = 4, & + URR_N_GAMMA = 5, & + URR_HEATING = 6 + ! Maximum number of partial fission reactions integer, parameter :: PARTIAL_FISSION_MAX = 4 diff --git a/src/cross_section.F90 b/src/cross_section.F90 index fde72e51b9..4332615abf 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -1107,6 +1107,25 @@ contains nuc % urr_data % multiply_smooth = .true. end if + ! if the inelastic competition flag indicates that the inelastic cross + ! section should be determined from a normal reaction cross section, we need + ! to set up a pointer to that reaction + nuc % urr_inelastic = NONE + if (nuc % urr_data % inelastic_flag > 0) then + do i = 1, nuc % n_reaction + if (nuc % reactions(i) % MT == nuc % urr_data % inelastic_flag) then + nuc % urr_inelastic = i + end if + end do + + ! Abort if no corresponding inelastic reaction was found + if (nuc % urr_inelastic == NONE) then + message = "Could not find inelastic reaction specified on " & + // "unresolved resonance probability table." + call fatal_error() + end if + end if + ! allocate incident energies and probability tables N = nuc % urr_data % n_energy M = nuc % urr_data % n_prob diff --git a/src/cross_section_header.F90 b/src/cross_section_header.F90 index 3eaf6d4462..2cc7b481c7 100644 --- a/src/cross_section_header.F90 +++ b/src/cross_section_header.F90 @@ -113,6 +113,7 @@ module cross_section_header ! Unresolved resonance data logical :: urr_present + integer :: urr_inelastic type(UrrData), pointer :: urr_data => null() ! Reactions diff --git a/src/global.F90 b/src/global.F90 index 4cfb6a6758..9c65efaf8a 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -73,6 +73,9 @@ module global integer :: n_grid ! number of points on unionized grid real(8), allocatable :: e_grid(:) ! energies on unionized grid + ! Unreoslved resonance probablity tables + logical :: urr_ptables_on = .false. + ! ============================================================================ ! TALLY-RELATED VARIABLES diff --git a/src/physics.F90 b/src/physics.F90 index 0db20f28c9..7f1487ddb6 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,7 +1,7 @@ module physics use constants - use cross_section_header, only: Nuclide, Reaction, DistEnergy + use cross_section_header, only: Nuclide, Reaction, DistEnergy, UrrData use endf, only: reaction_name, is_fission, is_scatter use error, only: fatal_error, warning use fission, only: nu_total, nu_prompt, nu_delayed @@ -279,6 +279,16 @@ contains if (index_sab > 0) call calculate_sab_xs(p, index_nuclide, index_sab) + ! if the particle is in the unresolved resonance range and there are + ! probability tables, we need to determine cross sections from the table + + if (urr_ptables_on .and. nuc % urr_present) then + if (p % E > nuc % urr_data % energy(1) .and. & + p % E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then + call calculate_urr_xs(p, index_nuclide) + end if + end if + ! Set last evaluated energy micro_xs(index_nuclide) % last_E = p % E @@ -368,6 +378,113 @@ contains end subroutine calculate_sab_xs +!=============================================================================== +! CALCULATE_URR_XS determines cross sections in the unresolved resonance range +! from probability tables +!=============================================================================== + + subroutine calculate_urr_xs(p, index_nuclide) + + type(Particle), pointer :: p + integer, intent(in) :: index_nuclide ! index into nuclides array + + integer :: i_energy ! index for energy + integer :: i_table ! index for table + real(8) :: f ! interpolation factor + real(8) :: r ! pseudo-random number + real(8) :: elastic ! smooth elastic cross section + real(8) :: absorption ! smooth absorption cross section + real(8) :: fission ! smooth fission cross section + real(8) :: inelastic ! smooth inelastic cross section + type(UrrData), pointer :: urr => null() + type(Nuclide), pointer :: nuc => null() + type(Reaction), pointer :: rxn => null() + + ! copy cross-sections already calculated + elastic = micro_xs(index_nuclide) % elastic + absorption = micro_xs(index_nuclide) % absorption + fission = micro_xs(index_nuclide) % fission + + ! get pointer to probability table + nuc => nuclides(index_nuclide) + urr => nuc % urr_data + + ! determine energy table + i_energy = 1 + do + if (p % E > urr % energy(i_energy)) exit + i_energy = i_energy + 1 + end do + + ! determine interpolation factor on table + f = (p % E - urr % energy(i_energy)) / & + (urr % energy(i_energy + 1) - urr % energy(i_energy)) + + ! sample probability table using the cumulative distribution + r = prn() + i_table = 1 + do + if (urr % prob(i_energy, URR_CUM_PROB, i_table) > r) exit + i_table = i_table + 1 + end do + + ! determine elastic, fission, and capture cross sections from probability + ! table + if (urr % interp == LINEAR_LINEAR) then + micro_xs(index_nuclide) % elastic = & + (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_table) + & + f * urr % prob(i_energy + 1, URR_ELASTIC, i_table) + micro_xs(index_nuclide) % fission = & + (ONE - f) * urr % prob(i_energy, URR_FISSION, i_table) + & + f * urr % prob(i_energy + 1, URR_FISSION, i_table) + micro_xs(index_nuclide) % absorption = & + micro_xs(index_nuclide) % fission + & + (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_table) + & + f * urr % prob(i_energy + 1, URR_N_GAMMA, i_table) + elseif (urr % interp == LOG_LOG) then + message = "Log-log interpolation on probability table not yet supported." + call fatal_error() + end if + + ! Determine treatment of inelastic scattering + inelastic = ZERO + if (urr % inelastic_flag > 0) then + ! Get pointer to inelastic scattering reaction + rxn => nuc % reactions(nuc % urr_inelastic) + + ! Get index on energy grid and interpolation factor + i_energy = micro_xs(index_nuclide) % index_grid + f = micro_xs(index_nuclide) % interp_factor + + ! Determine inelastic scattering cross section + if (i_energy >= rxn % IE) then + inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%IE + 1) + & + f * rxn % sigma(i_energy - rxn%IE + 2) + end if + end if + + ! Multiply by smooth cross-section if needed + if (urr % multiply_smooth) then + micro_xs(index_nuclide) % elastic = elastic * & + micro_xs(index_nuclide) % elastic + micro_xs(index_nuclide) % absorption = absorption * & + micro_xs(index_nuclide) % absorption + micro_xs(index_nuclide) % fission = fission * & + micro_xs(index_nuclide) % fission + end if + + ! Determine nu-fission cross section + if (nuc % fissionable) then + micro_xs(index_nuclide) % nu_fission = nu_total(nuc, p % E) * & + micro_xs(index_nuclide) % fission + end if + + ! Calculate total cross section as sum of partials + micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % elastic + & + inelastic + micro_xs(index_nuclide) % absorption + + end subroutine calculate_urr_xs + !=============================================================================== ! FIND_ENERGY_INDEX determines the index on the union energy grid and the ! interpolation factor for a particle at a certain energy From ffdcead9e4be5dfcfee52ffe0cdb67b8cf23b813 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 18 Dec 2011 17:59:29 -0500 Subject: [PATCH 20/33] Check for neighboring energy points that are the same on nuclide energy grid. --- src/physics.F90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/physics.F90 b/src/physics.F90 index 7f1487ddb6..f4dcbd28dd 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -235,8 +235,13 @@ contains ! TODO: If not using unionized energy grid, we need to find the index on the ! nuclide energy grid using lethargy mapping or whatever other technique - ! search nuclide energy grid + ! get index on nuclide energy grid IE = nuc % grid_index(p % IE) + + ! check for rare case where two energy points are the same + if (nuc % energy(IE) == nuc % energy(IE+1)) IE = IE + 1 + + ! calculate interpolation factor f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE)) micro_xs(index_nuclide) % index_grid = IE From 2898702569afa74e586e23f7c7d6e601fd250baa Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Dec 2011 09:56:43 -0500 Subject: [PATCH 21/33] Fixed cross section output. --- src/output.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/output.F90 b/src/output.F90 index d51f83ee09..8ebf16d846 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -780,7 +780,7 @@ contains size_energy = 0 end if - write(unit_,'(3X,A11,1X,F8.3,2X,I3,3X,I6,1X,I11,1X,I11)') & + write(unit_,'(3X,A11,1X,F8.3,2X,I4,2X,I6,1X,I11,1X,I11)') & reaction_name(rxn % MT), rxn % Q_value, rxn % TY, rxn % IE, & size_angle, size_energy From aed80ac08d6aa1f42249f93d52a2fdfbfea535b0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Dec 2011 12:50:08 -0500 Subject: [PATCH 22/33] Added check for unsupported external source type. --- src/source.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/source.F90 b/src/source.F90 index b3b9a552f1..dae40548a7 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -72,6 +72,10 @@ contains if (external_source%type == SRC_BOX) then p_min = external_source%values(1:3) p_max = external_source%values(4:6) + else + message = "Unsupported external source type: " // & + int_to_str(external_source%type) + call fatal_error() end if ! Initialize first cycle source bank From 7599ee7e311065df092ad38a74e9d4fda8205d62 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Dec 2011 13:10:47 -0500 Subject: [PATCH 23/33] Renamed cross_section modules to ace. --- src/DEPENDENCIES | 38 +++++++++---------- src/OBJECTS | 4 +- src/{cross_section.F90 => ace.F90} | 8 ++-- ...ross_section_header.F90 => ace_header.F90} | 4 +- src/fission.F90 | 10 ++--- src/global.F90 | 22 +++++------ src/initialize.F90 | 2 +- src/output.F90 | 16 ++++---- src/physics.F90 | 28 +++++++------- src/source.F90 | 18 ++++----- 10 files changed, 75 insertions(+), 75 deletions(-) rename src/{cross_section.F90 => ace.F90} (99%) rename src/{cross_section_header.F90 => ace_header.F90} (99%) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 8190c78b61..f385b0b940 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -1,18 +1,18 @@ -cross_section.o: constants.o -cross_section.o: cross_section_header.o -cross_section.o: datatypes.o -cross_section.o: datatypes_header.o -cross_section.o: endf.o -cross_section.o: error.o -cross_section.o: fileio.o -cross_section.o: fission.o -cross_section.o: global.o -cross_section.o: material_header.o -cross_section.o: output.o -cross_section.o: string.o +ace.o: ace_header.o +ace.o: constants.o +ace.o: datatypes.o +ace.o: datatypes_header.o +ace.o: endf.o +ace.o: error.o +ace.o: fileio.o +ace.o: fission.o +ace.o: global.o +ace.o: material_header.o +ace.o: output.o +ace.o: string.o -cross_section_header.o: constants.o -cross_section_header.o: endf_header.o +ace_header.o: constants.o +ace_header.o: endf_header.o datatypes.o: datatypes_header.o @@ -33,8 +33,8 @@ fileio.o: constants.o fileio.o: global.o fileio.o: string.o +fission.o: ace_header.o fission.o: constants.o -fission.o: cross_section_header.o fission.o: error.o fission.o: global.o fission.o: interpolation.o @@ -50,9 +50,9 @@ geometry.o: particle_header.o geometry.o: string.o geometry.o: tally.o +global.o: ace_header.o global.o: bank_header.o global.o: constants.o -global.o: cross_section_header.o global.o: datatypes_header.o global.o: geometry_header.o global.o: material_header.o @@ -62,8 +62,8 @@ global.o: source_header.o global.o: tally_header.o global.o: timing.o +initialize.o: ace.o initialize.o: constants.o -initialize.o: cross_section.o initialize.o: datatypes.o initialize.o: datatypes_header.o initialize.o: energy_grid.o @@ -145,8 +145,8 @@ output.o: tally_header.o particle_header.o: constants.o particle_header.o: geometry_header.o +physics.o: ace_header.o physics.o: constants.o -physics.o: cross_section_header.o physics.o: endf.o physics.o: error.o physics.o: fission.o @@ -172,9 +172,9 @@ search.o: constants.o search.o: error.o search.o: global.o +source.o: ace_header.o source.o: bank_header.o source.o: constants.o -source.o: cross_section_header.o source.o: error.o source.o: global.o source.o: output.o diff --git a/src/OBJECTS b/src/OBJECTS index 09f29ba6f8..6703d72439 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -1,7 +1,7 @@ objects = \ +ace.o \ +ace_header.o \ bank_header.o \ -cross_section_header.o \ -cross_section.o \ datatypes.o \ datatypes_header.o \ doppler.o \ diff --git a/src/cross_section.F90 b/src/ace.F90 similarity index 99% rename from src/cross_section.F90 rename to src/ace.F90 index 4332615abf..fe85df6aae 100644 --- a/src/cross_section.F90 +++ b/src/ace.F90 @@ -1,8 +1,8 @@ -module cross_section +module ace - use constants - use cross_section_header, only: Nuclide, Reaction, SAB_Table, XsListing, & + use ace_header, only: Nuclide, Reaction, SAB_Table, XsListing, & DistEnergy + use constants use datatypes, only: dict_create, dict_add_key, dict_get_key, & dict_has_key, dict_delete, dict_keys use datatypes_header, only: DictionaryCI, ListKeyValueCI @@ -1312,4 +1312,4 @@ contains end function get_real -end module cross_section +end module ace diff --git a/src/cross_section_header.F90 b/src/ace_header.F90 similarity index 99% rename from src/cross_section_header.F90 rename to src/ace_header.F90 index 2cc7b481c7..0c14a1ca5f 100644 --- a/src/cross_section_header.F90 +++ b/src/ace_header.F90 @@ -1,4 +1,4 @@ -module cross_section_header +module ace_header use constants, only: MAX_FILE_LEN use endf_header, only: Tab1 @@ -209,4 +209,4 @@ module cross_section_header real(8) :: nu_fission ! macroscopic production xs end type MaterialMacroXS -end module cross_section_header +end module ace_header diff --git a/src/fission.F90 b/src/fission.F90 index c544385361..f82bcf3540 100644 --- a/src/fission.F90 +++ b/src/fission.F90 @@ -1,11 +1,11 @@ module fission + use ace_header, only: Nuclide use constants - use cross_section_header, only: Nuclide - use error, only: fatal_error - use global, only: message - use interpolation, only: interpolate_tab1 - use search, only: binary_search + use error, only: fatal_error + use global, only: message + use interpolation, only: interpolate_tab1 + use search, only: binary_search implicit none diff --git a/src/global.F90 b/src/global.F90 index 9c65efaf8a..e904f64a29 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -1,17 +1,17 @@ module global - use bank_header, only: Bank + use ace_header, only: Nuclide, SAB_Table, xsListing, NuclideMicroXS, & + MaterialMacroXS + use bank_header, only: Bank use constants - use cross_section_header, only: Nuclide, SAB_Table, xsListing, & - NuclideMicroXS, MaterialMacroXS - use datatypes_header, only: DictionaryII, DictionaryCI - use geometry_header, only: Cell, Universe, Lattice, Surface - use material_header, only: Material - use mesh_header, only: StructuredMesh - use particle_header, only: Particle - use source_header, only: ExtSource - use tally_header, only: TallyObject, TallyMap - use timing, only: Timer + use datatypes_header, only: DictionaryII, DictionaryCI + use geometry_header, only: Cell, Universe, Lattice, Surface + use material_header, only: Material + use mesh_header, only: StructuredMesh + use particle_header, only: Particle + use source_header, only: ExtSource + use tally_header, only: TallyObject, TallyMap + use timing, only: Timer #ifdef MPI use mpi diff --git a/src/initialize.F90 b/src/initialize.F90 index 621df2c649..220c2a6809 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -1,7 +1,7 @@ module initialize + use ace, only: read_xs use constants - use cross_section, only: read_xs use datatypes, only: dict_create, dict_add_key, dict_get_key, & dict_has_key, dict_keys use datatypes_header, only: ListKeyValueII, DictionaryII diff --git a/src/output.F90 b/src/output.F90 index 8ebf16d846..32e964123c 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -2,16 +2,16 @@ module output use ISO_FORTRAN_ENV + use ace_header, only: Nuclide, Reaction, UrrData use constants - use cross_section_header, only: Nuclide, Reaction, UrrData - use datatypes, only: dict_get_key - use endf, only: reaction_name - use geometry_header, only: Cell, Universe, Surface + use datatypes, only: dict_get_key + use endf, only: reaction_name + use geometry_header, only: Cell, Universe, Surface use global - use mesh_header, only: StructuredMesh - use particle_header, only: Particle, LocalCoord - use string, only: upper_case, int_to_str, real_to_str - use tally_header, only: TallyObject + use mesh_header, only: StructuredMesh + use particle_header, only: Particle, LocalCoord + use string, only: upper_case, int_to_str, real_to_str + use tally_header, only: TallyObject implicit none diff --git a/src/physics.F90 b/src/physics.F90 index f4dcbd28dd..e9e86a862b 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,21 +1,21 @@ module physics + use ace_header, only: Nuclide, Reaction, DistEnergy, UrrData use constants - use cross_section_header, only: Nuclide, Reaction, DistEnergy, UrrData - use endf, only: reaction_name, is_fission, is_scatter - use error, only: fatal_error, warning - use fission, only: nu_total, nu_prompt, nu_delayed - use geometry, only: find_cell, distance_to_boundary, & - cross_surface, cross_lattice - use geometry_header, only: Universe, BASE_UNIVERSE + use endf, only: reaction_name, is_fission, is_scatter + use error, only: fatal_error, warning + use fission, only: nu_total, nu_prompt, 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 output, only: write_message - use particle_header, only: Particle, LocalCoord - use random_lcg, only: prn - use search, only: binary_search - use string, only: int_to_str, real_to_str - use tally, only: score_tally, score_surface_current + use interpolation, only: interpolate_tab1 + use output, only: write_message + use particle_header, only: Particle, LocalCoord + use random_lcg, only: prn + use search, only: binary_search + use string, only: int_to_str, real_to_str + use tally, only: score_tally, score_surface_current implicit none diff --git a/src/source.F90 b/src/source.F90 index dae40548a7..f296c2c5b6 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -1,15 +1,15 @@ module source - use bank_header, only: Bank - use constants, only: ONE, MAX_LINE_LEN - use cross_section_header, only: Nuclide - use error, only: fatal_error + use ace_header, only: Nuclide + use bank_header, only: Bank + use constants, only: ONE, MAX_LINE_LEN + use error, only: fatal_error use global - use output, only: write_message - use particle_header, only: Particle, initialize_particle - use physics, only: watt_spectrum - use random_lcg, only: prn, set_particle_seed - use string, only: int_to_str + use output, only: write_message + use particle_header, only: Particle, initialize_particle + use physics, only: watt_spectrum + use random_lcg, only: prn, set_particle_seed + use string, only: int_to_str implicit none From 1602fe672a8de40376bf09b22461139e00b083a0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Dec 2011 13:53:57 -0500 Subject: [PATCH 24/33] Moved cross section calculation routines into separate module. --- src/DEPENDENCIES | 11 ++ src/OBJECTS | 1 + src/cross_section.F90 | 408 ++++++++++++++++++++++++++++++++++++++++++ src/physics.F90 | 394 +--------------------------------------- 4 files changed, 422 insertions(+), 392 deletions(-) create mode 100644 src/cross_section.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index f385b0b940..977875c1cb 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -14,6 +14,16 @@ ace.o: string.o ace_header.o: constants.o ace_header.o: endf_header.o +cross_section.o: ace_header.o +cross_section.o: constants.o +cross_section.o: error.o +cross_section.o: fission.o +cross_section.o: global.o +cross_section.o: material_header.o +cross_section.o: particle_header.o +cross_section.o: random_lcg.o +cross_section.o: search.o + datatypes.o: datatypes_header.o doppler.o: constants.o @@ -147,6 +157,7 @@ particle_header.o: geometry_header.o physics.o: ace_header.o physics.o: constants.o +physics.o: cross_section.o physics.o: endf.o physics.o: error.o physics.o: fission.o diff --git a/src/OBJECTS b/src/OBJECTS index 6703d72439..bbf3b2a257 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -2,6 +2,7 @@ objects = \ ace.o \ ace_header.o \ bank_header.o \ +cross_section.o \ datatypes.o \ datatypes_header.o \ doppler.o \ diff --git a/src/cross_section.F90 b/src/cross_section.F90 new file mode 100644 index 0000000000..79ac18d00c --- /dev/null +++ b/src/cross_section.F90 @@ -0,0 +1,408 @@ +module cross_section + + use ace_header, only: Nuclide, SAB_Table, Reaction, UrrData + use constants + use error, only: fatal_error + use fission, only: nu_total + use global + use material_header, only: Material + use particle_header, only: Particle + use random_lcg, only: prn + use search, only: binary_search + + implicit none + +contains + +!=============================================================================== +! CALCULATE_XS determines the macroscopic cross sections for the material the +! particle is currently traveling through. +!=============================================================================== + + subroutine calculate_xs(p) + + type(Particle), pointer :: p + + integer :: i ! loop index over nuclides + integer :: index_nuclide ! index into nuclides array + integer :: index_sab ! index into sab_tables array + real(8) :: atom_density ! atom density of a nuclide + real(8) :: sab_threshold ! threshold for S(a,b) table + type(Material), pointer :: mat => null() ! current material + + ! Set all material macroscopic cross sections to zero + material_xs % total = ZERO + material_xs % elastic = ZERO + material_xs % absorption = ZERO + material_xs % fission = ZERO + material_xs % nu_fission = ZERO + + mat => materials(p % material) + + ! Find energy index on unionized grid + call find_energy_index(p) + + ! Check if there's an S(a,b) table for this material + if (mat % has_sab_table) then + sab_threshold = sab_tables(mat % sab_table) % threshold_inelastic + else + sab_threshold = ZERO + end if + + ! Add contribution from each nuclide in material + do i = 1, mat % n_nuclides + ! Determine microscopic cross sections for this nuclide + index_nuclide = mat % nuclide(i) + + ! Determine whether to use S(a,b) based on energy of particle and whether + ! the nuclide matches the S(a,b) table + if (p % E < sab_threshold .and. i == mat % sab_nuclide) then + index_sab = mat % sab_table + else + index_sab = 0 + end if + + ! Calculate microscopic cross section for this nuclide + if (p % E /= micro_xs(index_nuclide) % last_E) then + call calculate_nuclide_xs(p, index_nuclide, index_sab) + end if + + ! Copy atom density of nuclide in material + atom_density = mat % atom_density(i) + + ! Add contributions to material macroscopic total cross section + material_xs % total = material_xs % total + & + atom_density * micro_xs(index_nuclide) % total + + ! Add contributions to material macroscopic scattering cross section + material_xs % elastic = material_xs % elastic + & + atom_density * micro_xs(index_nuclide) % elastic + + ! Add contributions to material macroscopic absorption cross section + material_xs % absorption = material_xs % absorption + & + atom_density * micro_xs(index_nuclide) % absorption + + ! Add contributions to material macroscopic fission cross section + material_xs % fission = material_xs % fission + & + atom_density * micro_xs(index_nuclide) % fission + + ! Add contributions to material macroscopic nu-fission cross section + material_xs % nu_fission = material_xs % nu_fission + & + atom_density * micro_xs(index_nuclide) % nu_fission + end do + + end subroutine calculate_xs + +!=============================================================================== +! CALCULATE_NUCLIDE_XS determines microscopic cross sections for a nuclide of a +! given index in the nuclides array at the energy of the given particle +!=============================================================================== + + subroutine calculate_nuclide_xs(p, index_nuclide, index_sab) + + type(Particle), pointer :: p + integer, intent(in) :: index_nuclide ! index into nuclides array + integer, intent(in) :: index_sab ! index into sab_tables array + + integer :: IE ! index on nuclide energy grid + real(8) :: f ! interp factor on nuclide energy grid + type(Nuclide), pointer :: nuc => null() + + ! Set pointer to nuclide + nuc => nuclides(index_nuclide) + + ! TODO: If not using unionized energy grid, we need to find the index on the + ! nuclide energy grid using lethargy mapping or whatever other technique + + ! get index on nuclide energy grid + IE = nuc % grid_index(p % IE) + + ! check for rare case where two energy points are the same + if (nuc % energy(IE) == nuc % energy(IE+1)) IE = IE + 1 + + ! calculate interpolation factor + f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE)) + + micro_xs(index_nuclide) % index_grid = IE + micro_xs(index_nuclide) % interp_factor = f + + ! Initialize sab treatment to false + micro_xs(index_nuclide) % use_sab = .false. + micro_xs(index_nuclide) % elastic_sab = ZERO + + ! Initialize nuclide cross-sections to zero + micro_xs(index_nuclide) % fission = ZERO + micro_xs(index_nuclide) % nu_fission = ZERO + + ! Calculate microscopic nuclide total cross section + micro_xs(index_nuclide) % total = & + (ONE-f) * nuc % total(IE) + f * nuc % total(IE+1) + + ! Calculate microscopic nuclide total cross section + micro_xs(index_nuclide) % elastic = & + (ONE-f) * nuc % elastic(IE) + f * nuc % elastic(IE+1) + + ! Calculate microscopic nuclide absorption cross section + micro_xs(index_nuclide) % absorption = & + (ONE-f) * nuc % absorption(IE) + f * nuc % absorption(IE+1) + + if (nuc % fissionable) then + ! Calculate microscopic nuclide total cross section + micro_xs(index_nuclide) % fission = & + (ONE-f) * nuc % fission(IE) + f * nuc % fission(IE+1) + + ! Calculate microscopic nuclide nu-fission cross section + micro_xs(index_nuclide) % nu_fission = & + (ONE-f) * nuc % nu_fission(IE) + f * nuc % nu_fission(IE+1) + end if + + ! If there is S(a,b) data for this nuclide, we need to do a few + ! things. Since the total cross section was based on non-S(a,b) data, we + ! need to correct it by subtracting the non-S(a,b) elastic cross section and + ! then add back in the calculated S(a,b) elastic+inelastic cross section. + + if (index_sab > 0) call calculate_sab_xs(p, index_nuclide, index_sab) + + ! if the particle is in the unresolved resonance range and there are + ! probability tables, we need to determine cross sections from the table + + if (urr_ptables_on .and. nuc % urr_present) then + if (p % E > nuc % urr_data % energy(1) .and. & + p % E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then + call calculate_urr_xs(p, index_nuclide) + end if + end if + + ! Set last evaluated energy + micro_xs(index_nuclide) % last_E = p % E + + end subroutine calculate_nuclide_xs + +!=============================================================================== +! CALCULATE_SAB_XS determines the elastic and inelastic scattering +! cross-sections in the thermal energy range. These cross sections replace +! whatever data were taken from the normal Nuclide table. +!=============================================================================== + + subroutine calculate_sab_xs(p, index_nuclide, index_sab) + + type(Particle), pointer :: p + integer, intent(in) :: index_nuclide ! index into nuclides array + integer, intent(in) :: index_sab ! index into sab_tables array + + integer :: IE_sab ! index on S(a,b) energy grid + real(8) :: f_sab ! interp factor on S(a,b) energy grid + real(8) :: inelastic ! S(a,b) inelastic cross section + real(8) :: elastic ! S(a,b) elastic cross section + type(SAB_Table), pointer :: sab => null() + + ! Set flag that S(a,b) treatment should be used for scattering + micro_xs(index_nuclide) % use_sab = .true. + + ! Get pointer to S(a,b) table + sab => sab_tables(index_sab) + + ! Get index and interpolation factor for inelastic grid + if (p%E < sab % inelastic_e_in(1)) then + IE_sab = 1 + f_sab = ZERO + else + IE_sab = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E) + f_sab = (p%E - sab%inelastic_e_in(IE_sab)) / & + (sab%inelastic_e_in(IE_sab+1) - sab%inelastic_e_in(IE_sab)) + end if + + ! Calculate S(a,b) inelastic scattering cross section + inelastic = (ONE-f_sab) * sab % inelastic_sigma(IE_sab) + f_sab * & + sab % inelastic_sigma(IE_sab + 1) + + ! Check for elastic data + if (p % E < sab % threshold_elastic) then + ! Determine whether elastic scattering is given in the coherent or + ! incoherent approximation. For coherent, the cross section is + ! represented as P/E whereas for incoherent, it is simply P + + if (sab % elastic_mode == SAB_ELASTIC_EXACT) then + if (p % E < sab % elastic_e_in(1)) then + ! If energy is below that of the lowest Bragg peak, the elastic + ! cross section will be zero + elastic = ZERO + else + IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) + elastic = sab % elastic_P(IE_sab) / p % E + end if + else + ! Determine index on elastic energy grid + if (p % E < sab % elastic_e_in(1)) then + IE_sab = 1 + else + IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) + end if + + ! Get interpolation factor for elastic grid + f_sab = (p%E - sab%elastic_e_in(IE_sab))/(sab%elastic_e_in(IE_sab+1) - & + sab%elastic_e_in(IE_sab)) + + ! Calculate S(a,b) elastic scattering cross section + elastic = (ONE-f_sab) * sab % elastic_P(IE_sab) + f_sab * & + sab % elastic_P(IE_sab + 1) + end if + else + ! No elastic data + elastic = ZERO + end if + + ! Correct total and elastic cross sections + micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % total - & + micro_xs(index_nuclide) % elastic + inelastic + elastic + micro_xs(index_nuclide) % elastic = inelastic + elastic + + ! Store S(a,b) elastic cross section for sampling later + micro_xs(index_nuclide) % elastic_sab = elastic + + end subroutine calculate_sab_xs + +!=============================================================================== +! CALCULATE_URR_XS determines cross sections in the unresolved resonance range +! from probability tables +!=============================================================================== + + subroutine calculate_urr_xs(p, index_nuclide) + + type(Particle), pointer :: p + integer, intent(in) :: index_nuclide ! index into nuclides array + + integer :: i_energy ! index for energy + integer :: i_table ! index for table + real(8) :: f ! interpolation factor + real(8) :: r ! pseudo-random number + real(8) :: elastic ! smooth elastic cross section + real(8) :: absorption ! smooth absorption cross section + real(8) :: fission ! smooth fission cross section + real(8) :: inelastic ! smooth inelastic cross section + type(UrrData), pointer :: urr => null() + type(Nuclide), pointer :: nuc => null() + type(Reaction), pointer :: rxn => null() + + ! copy cross-sections already calculated + elastic = micro_xs(index_nuclide) % elastic + absorption = micro_xs(index_nuclide) % absorption + fission = micro_xs(index_nuclide) % fission + + ! get pointer to probability table + nuc => nuclides(index_nuclide) + urr => nuc % urr_data + + ! determine energy table + i_energy = 1 + do + if (p % E > urr % energy(i_energy)) exit + i_energy = i_energy + 1 + end do + + ! determine interpolation factor on table + f = (p % E - urr % energy(i_energy)) / & + (urr % energy(i_energy + 1) - urr % energy(i_energy)) + + ! sample probability table using the cumulative distribution + r = prn() + i_table = 1 + do + if (urr % prob(i_energy, URR_CUM_PROB, i_table) > r) exit + i_table = i_table + 1 + end do + + ! determine elastic, fission, and capture cross sections from probability + ! table + if (urr % interp == LINEAR_LINEAR) then + micro_xs(index_nuclide) % elastic = & + (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_table) + & + f * urr % prob(i_energy + 1, URR_ELASTIC, i_table) + micro_xs(index_nuclide) % fission = & + (ONE - f) * urr % prob(i_energy, URR_FISSION, i_table) + & + f * urr % prob(i_energy + 1, URR_FISSION, i_table) + micro_xs(index_nuclide) % absorption = & + micro_xs(index_nuclide) % fission + & + (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_table) + & + f * urr % prob(i_energy + 1, URR_N_GAMMA, i_table) + elseif (urr % interp == LOG_LOG) then + message = "Log-log interpolation on probability table not yet supported." + call fatal_error() + end if + + ! Determine treatment of inelastic scattering + inelastic = ZERO + if (urr % inelastic_flag > 0) then + ! Get pointer to inelastic scattering reaction + rxn => nuc % reactions(nuc % urr_inelastic) + + ! Get index on energy grid and interpolation factor + i_energy = micro_xs(index_nuclide) % index_grid + f = micro_xs(index_nuclide) % interp_factor + + ! Determine inelastic scattering cross section + if (i_energy >= rxn % IE) then + inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%IE + 1) + & + f * rxn % sigma(i_energy - rxn%IE + 2) + end if + end if + + ! Multiply by smooth cross-section if needed + if (urr % multiply_smooth) then + micro_xs(index_nuclide) % elastic = elastic * & + micro_xs(index_nuclide) % elastic + micro_xs(index_nuclide) % absorption = absorption * & + micro_xs(index_nuclide) % absorption + micro_xs(index_nuclide) % fission = fission * & + micro_xs(index_nuclide) % fission + end if + + ! Determine nu-fission cross section + if (nuc % fissionable) then + micro_xs(index_nuclide) % nu_fission = nu_total(nuc, p % E) * & + micro_xs(index_nuclide) % fission + end if + + ! Calculate total cross section as sum of partials + micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % elastic + & + inelastic + micro_xs(index_nuclide) % absorption + + end subroutine calculate_urr_xs + +!=============================================================================== +! FIND_ENERGY_INDEX determines the index on the union energy grid and the +! interpolation factor for a particle at a certain energy +!=============================================================================== + + subroutine find_energy_index(p) + + type(Particle), pointer :: p + + integer :: IE ! index on union energy grid + real(8) :: E ! energy of particle + real(8) :: interp ! interpolation factor + + ! copy particle's energy + E = p % E + + ! if particle's energy is outside of energy grid range, set to first or last + ! index. Otherwise, do a binary search through the union energy grid. + if (E < e_grid(1)) then + IE = 1 + elseif (E > e_grid(n_grid)) then + IE = n_grid - 1 + else + IE = binary_search(e_grid, n_grid, E) + end if + + ! calculate the interpolation factor -- note this will be outside of [0,1) + ! for a particle outside the energy range of the union grid + interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE)) + + ! set particle attributes + p % IE = IE + p % interp = interp + + end subroutine find_energy_index + +end module cross_section diff --git a/src/physics.F90 b/src/physics.F90 index e9e86a862b..86650d3141 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,7 +1,8 @@ module physics - use ace_header, only: Nuclide, Reaction, DistEnergy, UrrData + use ace_header, only: Nuclide, Reaction, DistEnergy use constants + use cross_section, only: calculate_xs, find_energy_index use endf, only: reaction_name, is_fission, is_scatter use error, only: fatal_error, warning use fission, only: nu_total, nu_prompt, nu_delayed @@ -135,397 +136,6 @@ contains end subroutine transport -!=============================================================================== -! CALCULATE_XS determines the macroscopic cross sections for the material the -! particle is currently traveling through. -!=============================================================================== - - subroutine calculate_xs(p) - - type(Particle), pointer :: p - - integer :: i ! loop index over nuclides - integer :: index_nuclide ! index into nuclides array - integer :: index_sab ! index into sab_tables array - real(8) :: atom_density ! atom density of a nuclide - real(8) :: sab_threshold ! threshold for S(a,b) table - type(Material), pointer :: mat => null() ! current material - - ! Set all material macroscopic cross sections to zero - material_xs % total = ZERO - material_xs % elastic = ZERO - material_xs % absorption = ZERO - material_xs % fission = ZERO - material_xs % nu_fission = ZERO - - mat => materials(p % material) - - ! Find energy index on unionized grid - call find_energy_index(p) - - ! Check if there's an S(a,b) table for this material - if (mat % has_sab_table) then - sab_threshold = sab_tables(mat % sab_table) % threshold_inelastic - else - sab_threshold = ZERO - end if - - ! Add contribution from each nuclide in material - do i = 1, mat % n_nuclides - ! Determine microscopic cross sections for this nuclide - index_nuclide = mat % nuclide(i) - - ! Determine whether to use S(a,b) based on energy of particle and whether - ! the nuclide matches the S(a,b) table - if (p % E < sab_threshold .and. i == mat % sab_nuclide) then - index_sab = mat % sab_table - else - index_sab = 0 - end if - - ! Calculate microscopic cross section for this nuclide - if (p % E /= micro_xs(index_nuclide) % last_E) then - call calculate_nuclide_xs(p, index_nuclide, index_sab) - end if - - ! Copy atom density of nuclide in material - atom_density = mat % atom_density(i) - - ! Add contributions to material macroscopic total cross section - material_xs % total = material_xs % total + & - atom_density * micro_xs(index_nuclide) % total - - ! Add contributions to material macroscopic scattering cross section - material_xs % elastic = material_xs % elastic + & - atom_density * micro_xs(index_nuclide) % elastic - - ! Add contributions to material macroscopic absorption cross section - material_xs % absorption = material_xs % absorption + & - atom_density * micro_xs(index_nuclide) % absorption - - ! Add contributions to material macroscopic fission cross section - material_xs % fission = material_xs % fission + & - atom_density * micro_xs(index_nuclide) % fission - - ! Add contributions to material macroscopic nu-fission cross section - material_xs % nu_fission = material_xs % nu_fission + & - atom_density * micro_xs(index_nuclide) % nu_fission - end do - - end subroutine calculate_xs - -!=============================================================================== -! CALCULATE_NUCLIDE_XS determines microscopic cross sections for a nuclide of a -! given index in the nuclides array at the energy of the given particle -!=============================================================================== - - subroutine calculate_nuclide_xs(p, index_nuclide, index_sab) - - type(Particle), pointer :: p - integer, intent(in) :: index_nuclide ! index into nuclides array - integer, intent(in) :: index_sab ! index into sab_tables array - - integer :: IE ! index on nuclide energy grid - real(8) :: f ! interp factor on nuclide energy grid - type(Nuclide), pointer :: nuc => null() - - ! Set pointer to nuclide - nuc => nuclides(index_nuclide) - - ! TODO: If not using unionized energy grid, we need to find the index on the - ! nuclide energy grid using lethargy mapping or whatever other technique - - ! get index on nuclide energy grid - IE = nuc % grid_index(p % IE) - - ! check for rare case where two energy points are the same - if (nuc % energy(IE) == nuc % energy(IE+1)) IE = IE + 1 - - ! calculate interpolation factor - f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE)) - - micro_xs(index_nuclide) % index_grid = IE - micro_xs(index_nuclide) % interp_factor = f - - ! Initialize sab treatment to false - micro_xs(index_nuclide) % use_sab = .false. - micro_xs(index_nuclide) % elastic_sab = ZERO - - ! Initialize nuclide cross-sections to zero - micro_xs(index_nuclide) % fission = ZERO - micro_xs(index_nuclide) % nu_fission = ZERO - - ! Calculate microscopic nuclide total cross section - micro_xs(index_nuclide) % total = & - (ONE-f) * nuc % total(IE) + f * nuc % total(IE+1) - - ! Calculate microscopic nuclide total cross section - micro_xs(index_nuclide) % elastic = & - (ONE-f) * nuc % elastic(IE) + f * nuc % elastic(IE+1) - - ! Calculate microscopic nuclide absorption cross section - micro_xs(index_nuclide) % absorption = & - (ONE-f) * nuc % absorption(IE) + f * nuc % absorption(IE+1) - - if (nuc % fissionable) then - ! Calculate microscopic nuclide total cross section - micro_xs(index_nuclide) % fission = & - (ONE-f) * nuc % fission(IE) + f * nuc % fission(IE+1) - - ! Calculate microscopic nuclide nu-fission cross section - micro_xs(index_nuclide) % nu_fission = & - (ONE-f) * nuc % nu_fission(IE) + f * nuc % nu_fission(IE+1) - end if - - ! If there is S(a,b) data for this nuclide, we need to do a few - ! things. Since the total cross section was based on non-S(a,b) data, we - ! need to correct it by subtracting the non-S(a,b) elastic cross section and - ! then add back in the calculated S(a,b) elastic+inelastic cross section. - - if (index_sab > 0) call calculate_sab_xs(p, index_nuclide, index_sab) - - ! if the particle is in the unresolved resonance range and there are - ! probability tables, we need to determine cross sections from the table - - if (urr_ptables_on .and. nuc % urr_present) then - if (p % E > nuc % urr_data % energy(1) .and. & - p % E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then - call calculate_urr_xs(p, index_nuclide) - end if - end if - - ! Set last evaluated energy - micro_xs(index_nuclide) % last_E = p % E - - end subroutine calculate_nuclide_xs - -!=============================================================================== -! CALCULATE_SAB_XS determines the elastic and inelastic scattering -! cross-sections in the thermal energy range. These cross sections replace -! whatever data were taken from the normal Nuclide table. -!=============================================================================== - - subroutine calculate_sab_xs(p, index_nuclide, index_sab) - - type(Particle), pointer :: p - integer, intent(in) :: index_nuclide ! index into nuclides array - integer, intent(in) :: index_sab ! index into sab_tables array - - integer :: IE_sab ! index on S(a,b) energy grid - real(8) :: f_sab ! interp factor on S(a,b) energy grid - real(8) :: inelastic ! S(a,b) inelastic cross section - real(8) :: elastic ! S(a,b) elastic cross section - type(SAB_Table), pointer :: sab => null() - - ! Set flag that S(a,b) treatment should be used for scattering - micro_xs(index_nuclide) % use_sab = .true. - - ! Get pointer to S(a,b) table - sab => sab_tables(index_sab) - - ! Get index and interpolation factor for inelastic grid - if (p%E < sab % inelastic_e_in(1)) then - IE_sab = 1 - f_sab = ZERO - else - IE_sab = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E) - f_sab = (p%E - sab%inelastic_e_in(IE_sab)) / & - (sab%inelastic_e_in(IE_sab+1) - sab%inelastic_e_in(IE_sab)) - end if - - ! Calculate S(a,b) inelastic scattering cross section - inelastic = (ONE-f_sab) * sab % inelastic_sigma(IE_sab) + f_sab * & - sab % inelastic_sigma(IE_sab + 1) - - ! Check for elastic data - if (p % E < sab % threshold_elastic) then - ! Determine whether elastic scattering is given in the coherent or - ! incoherent approximation. For coherent, the cross section is - ! represented as P/E whereas for incoherent, it is simply P - - if (sab % elastic_mode == SAB_ELASTIC_EXACT) then - if (p % E < sab % elastic_e_in(1)) then - ! If energy is below that of the lowest Bragg peak, the elastic - ! cross section will be zero - elastic = ZERO - else - IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) - elastic = sab % elastic_P(IE_sab) / p % E - end if - else - ! Determine index on elastic energy grid - if (p % E < sab % elastic_e_in(1)) then - IE_sab = 1 - else - IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) - end if - - ! Get interpolation factor for elastic grid - f_sab = (p%E - sab%elastic_e_in(IE_sab))/(sab%elastic_e_in(IE_sab+1) - & - sab%elastic_e_in(IE_sab)) - - ! Calculate S(a,b) elastic scattering cross section - elastic = (ONE-f_sab) * sab % elastic_P(IE_sab) + f_sab * & - sab % elastic_P(IE_sab + 1) - end if - else - ! No elastic data - elastic = ZERO - end if - - ! Correct total and elastic cross sections - micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % total - & - micro_xs(index_nuclide) % elastic + inelastic + elastic - micro_xs(index_nuclide) % elastic = inelastic + elastic - - ! Store S(a,b) elastic cross section for sampling later - micro_xs(index_nuclide) % elastic_sab = elastic - - end subroutine calculate_sab_xs - -!=============================================================================== -! CALCULATE_URR_XS determines cross sections in the unresolved resonance range -! from probability tables -!=============================================================================== - - subroutine calculate_urr_xs(p, index_nuclide) - - type(Particle), pointer :: p - integer, intent(in) :: index_nuclide ! index into nuclides array - - integer :: i_energy ! index for energy - integer :: i_table ! index for table - real(8) :: f ! interpolation factor - real(8) :: r ! pseudo-random number - real(8) :: elastic ! smooth elastic cross section - real(8) :: absorption ! smooth absorption cross section - real(8) :: fission ! smooth fission cross section - real(8) :: inelastic ! smooth inelastic cross section - type(UrrData), pointer :: urr => null() - type(Nuclide), pointer :: nuc => null() - type(Reaction), pointer :: rxn => null() - - ! copy cross-sections already calculated - elastic = micro_xs(index_nuclide) % elastic - absorption = micro_xs(index_nuclide) % absorption - fission = micro_xs(index_nuclide) % fission - - ! get pointer to probability table - nuc => nuclides(index_nuclide) - urr => nuc % urr_data - - ! determine energy table - i_energy = 1 - do - if (p % E > urr % energy(i_energy)) exit - i_energy = i_energy + 1 - end do - - ! determine interpolation factor on table - f = (p % E - urr % energy(i_energy)) / & - (urr % energy(i_energy + 1) - urr % energy(i_energy)) - - ! sample probability table using the cumulative distribution - r = prn() - i_table = 1 - do - if (urr % prob(i_energy, URR_CUM_PROB, i_table) > r) exit - i_table = i_table + 1 - end do - - ! determine elastic, fission, and capture cross sections from probability - ! table - if (urr % interp == LINEAR_LINEAR) then - micro_xs(index_nuclide) % elastic = & - (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_table) + & - f * urr % prob(i_energy + 1, URR_ELASTIC, i_table) - micro_xs(index_nuclide) % fission = & - (ONE - f) * urr % prob(i_energy, URR_FISSION, i_table) + & - f * urr % prob(i_energy + 1, URR_FISSION, i_table) - micro_xs(index_nuclide) % absorption = & - micro_xs(index_nuclide) % fission + & - (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_table) + & - f * urr % prob(i_energy + 1, URR_N_GAMMA, i_table) - elseif (urr % interp == LOG_LOG) then - message = "Log-log interpolation on probability table not yet supported." - call fatal_error() - end if - - ! Determine treatment of inelastic scattering - inelastic = ZERO - if (urr % inelastic_flag > 0) then - ! Get pointer to inelastic scattering reaction - rxn => nuc % reactions(nuc % urr_inelastic) - - ! Get index on energy grid and interpolation factor - i_energy = micro_xs(index_nuclide) % index_grid - f = micro_xs(index_nuclide) % interp_factor - - ! Determine inelastic scattering cross section - if (i_energy >= rxn % IE) then - inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%IE + 1) + & - f * rxn % sigma(i_energy - rxn%IE + 2) - end if - end if - - ! Multiply by smooth cross-section if needed - if (urr % multiply_smooth) then - micro_xs(index_nuclide) % elastic = elastic * & - micro_xs(index_nuclide) % elastic - micro_xs(index_nuclide) % absorption = absorption * & - micro_xs(index_nuclide) % absorption - micro_xs(index_nuclide) % fission = fission * & - micro_xs(index_nuclide) % fission - end if - - ! Determine nu-fission cross section - if (nuc % fissionable) then - micro_xs(index_nuclide) % nu_fission = nu_total(nuc, p % E) * & - micro_xs(index_nuclide) % fission - end if - - ! Calculate total cross section as sum of partials - micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % elastic + & - inelastic + micro_xs(index_nuclide) % absorption - - end subroutine calculate_urr_xs - -!=============================================================================== -! FIND_ENERGY_INDEX determines the index on the union energy grid and the -! interpolation factor for a particle at a certain energy -!=============================================================================== - - subroutine find_energy_index(p) - - type(Particle), pointer :: p - - integer :: IE ! index on union energy grid - real(8) :: E ! energy of particle - real(8) :: interp ! interpolation factor - - ! copy particle's energy - E = p % E - - ! if particle's energy is outside of energy grid range, set to first or last - ! index. Otherwise, do a binary search through the union energy grid. - if (E < e_grid(1)) then - IE = 1 - elseif (E > e_grid(n_grid)) then - IE = n_grid - 1 - else - IE = binary_search(e_grid, n_grid, E) - end if - - ! calculate the interpolation factor -- note this will be outside of [0,1) - ! for a particle outside the energy range of the union grid - interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE)) - - ! set particle attributes - p % IE = IE - p % interp = interp - - end subroutine find_energy_index - !=============================================================================== ! COLLISION samples a nuclide and reaction and then calls the appropriate ! routine for that reaction From a9fc06987e08f91421de0a0560fcb387da1c9b29 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Dec 2011 13:56:16 -0500 Subject: [PATCH 25/33] Removed redundant call to find_energy_index. --- src/physics.F90 | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 86650d3141..2a552c875b 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -2,7 +2,7 @@ module physics use ace_header, only: Nuclide, Reaction, DistEnergy use constants - use cross_section, only: calculate_xs, find_energy_index + use cross_section, only: calculate_xs use endf, only: reaction_name, is_fission, is_scatter use error, only: fatal_error, warning use fission, only: nu_total, nu_prompt, nu_delayed @@ -198,9 +198,6 @@ contains ! Reset number of particles banked during collision p % n_bank = 0 - ! find energy index, interpolation factor - call find_energy_index(p) - end subroutine collision !=============================================================================== From 154da86236ec801a8fe975ada23d3ce7127f0c2e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Dec 2011 14:23:02 -0500 Subject: [PATCH 26/33] Changed many line continuations. --- src/ace.F90 | 6 +-- src/constants.F90 | 96 ++++++++++++++++++++--------------------- src/geometry.F90 | 8 ++-- src/geometry_header.F90 | 4 +- src/initialize.F90 | 20 ++++----- src/input_xml.F90 | 48 ++++++++++----------- src/intercycle.F90 | 2 +- src/mpi_routines.F90 | 24 +++++------ src/output.F90 | 4 +- src/physics.F90 | 24 +++++------ src/string.F90 | 2 +- 11 files changed, 118 insertions(+), 120 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index fe85df6aae..f4d5cf374e 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -203,8 +203,8 @@ contains ! Check to make sure S(a,b) table matched a nuclide if (mat % sab_nuclide == 0) then - message = "S(a,b) table " // trim(mat % sab_name) // " did not match " & - // "any nuclide on material " // trim(int_to_str(mat % id)) + message = "S(a,b) table " // trim(mat % sab_name) // " did not match & + &any nuclide on material " // trim(int_to_str(mat % id)) call fatal_error() end if end if @@ -287,7 +287,7 @@ contains read(UNIT=in, FMT=100) comment, mat, & (zaids(i), awrs(i), i=1,16), NXS, JXS 100 format(A70,A10/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/& - &,8I9/8I9/8I9/8I9/8I9/8I9) + ,8I9/8I9/8I9/8I9/8I9/8I9) ! determine table length length = NXS(1) diff --git a/src/constants.F90 b/src/constants.F90 index 5188035f83..1f58ca4918 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -54,64 +54,64 @@ module constants ! GEOMETRY-RELATED CONSTANTS ! Boundary conditions - integer, parameter :: & - & BC_TRANSMIT = 0, & ! Transmission boundary condition (default) - & BC_VACUUM = 1, & ! Vacuum boundary condition - & BC_REFLECT = 2, & ! Reflecting boundary condition - & BC_PERIODIC = 3 ! Periodic boundary condition + integer, parameter :: & + BC_TRANSMIT = 0, & ! Transmission boundary condition (default) + BC_VACUUM = 1, & ! Vacuum boundary condition + BC_REFLECT = 2, & ! Reflecting boundary condition + BC_PERIODIC = 3 ! Periodic boundary condition ! Logical operators for cell definitions - integer, parameter :: & - & OP_LEFT_PAREN = huge(0), & ! Left parentheses - & OP_RIGHT_PAREN = huge(0) - 1, & ! Right parentheses - & OP_UNION = huge(0) - 2, & ! Union operator - & OP_DIFFERENCE = huge(0) - 3 ! Difference operator + integer, parameter :: & + OP_LEFT_PAREN = huge(0), & ! Left parentheses + OP_RIGHT_PAREN = huge(0) - 1, & ! Right parentheses + OP_UNION = huge(0) - 2, & ! Union operator + OP_DIFFERENCE = huge(0) - 3 ! Difference operator ! Cell types - integer, parameter :: & - & CELL_NORMAL = 1, & ! Cell with a specified material - & CELL_FILL = 2, & ! Cell filled by a separate universe - & CELL_LATTICE = 3, & ! Cell filled with a lattice - & CELL_VOID = -1 + integer, parameter :: & + CELL_NORMAL = 1, & ! Cell with a specified material + CELL_FILL = 2, & ! Cell filled by a separate universe + CELL_LATTICE = 3, & ! Cell filled with a lattice + CELL_VOID = -1 ! Lattice types - integer, parameter :: & - & LATTICE_RECT = 1, & - & LATTICE_HEX = 2 + integer, parameter :: & + LATTICE_RECT = 1, & + LATTICE_HEX = 2 ! Lattice boundary crossings - integer, parameter :: & - & LATTICE_LEFT = 1, & - & LATTICE_RIGHT = 2, & - & LATTICE_BOTTOM = 3, & - & LATTICE_TOP = 4 + integer, parameter :: & + LATTICE_LEFT = 1, & + LATTICE_RIGHT = 2, & + LATTICE_BOTTOM = 3, & + LATTICE_TOP = 4 ! Surface types - integer, parameter :: & - & SURF_PX = 1, & ! Plane parallel to x-plane - & SURF_PY = 2, & ! Plane parallel to y-plane - & SURF_PZ = 3, & ! Plane parallel to z-plane - & SURF_PLANE = 4, & ! Arbitrary plane - & SURF_CYL_X = 5, & ! Cylinder along x-axis - & SURF_CYL_Y = 6, & ! Cylinder along y-axis - & SURF_CYL_Z = 7, & ! Cylinder along z-axis - & SURF_SPHERE = 8, & ! Sphere - & SURF_BOX_X = 9, & ! Box extending infinitely in x-direction - & SURF_BOX_Y = 10, & ! Box extending infinitely in y-direction - & SURF_BOX_Z = 11, & ! Box extending infinitely in z-direction - & SURF_BOX = 12, & ! Rectangular prism - & SURF_GQ = 13 ! General quadratic surface + integer, parameter :: & + SURF_PX = 1, & ! Plane parallel to x-plane + SURF_PY = 2, & ! Plane parallel to y-plane + SURF_PZ = 3, & ! Plane parallel to z-plane + SURF_PLANE = 4, & ! Arbitrary plane + SURF_CYL_X = 5, & ! Cylinder along x-axis + SURF_CYL_Y = 6, & ! Cylinder along y-axis + SURF_CYL_Z = 7, & ! Cylinder along z-axis + SURF_SPHERE = 8, & ! Sphere + SURF_BOX_X = 9, & ! Box extending infinitely in x-direction + SURF_BOX_Y = 10, & ! Box extending infinitely in y-direction + SURF_BOX_Z = 11, & ! Box extending infinitely in z-direction + SURF_BOX = 12, & ! Rectangular prism + SURF_GQ = 13 ! General quadratic surface ! Surface senses - integer, parameter :: & - & SENSE_POSITIVE = 1, & - & SENSE_NEGATIVE = -1 + integer, parameter :: & + SENSE_POSITIVE = 1, & + SENSE_NEGATIVE = -1 ! ============================================================================ ! CROSS SECTION RELATED CONSTANTS ! Interpolation flag - integer, parameter :: & + integer, parameter :: & HISTOGRAM = 1, & ! y is constant in x LINEAR_LINEAR = 2, & ! y is linear in x LINEAR_LOG = 3, & ! y is linear in ln(x) @@ -126,19 +126,19 @@ module constants ! Angular distribution type integer, parameter :: & - & ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution - & ANGLE_32_EQUI = 2, & ! 32 equiprobable bins - & ANGLE_TABULAR = 3 ! Tabular angular distribution + ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution + ANGLE_32_EQUI = 2, & ! 32 equiprobable bins + ANGLE_TABULAR = 3 ! Tabular angular distribution ! Secondary energy mode for S(a,b) inelastic scattering integer, parameter :: & - & SAB_SECONDARY_EQUAL = 0, & ! Equally-likely outgoing energy bins - & SAB_SECONDARY_SKEWED = 1 ! Skewed outgoing energy bins + SAB_SECONDARY_EQUAL = 0, & ! Equally-likely outgoing energy bins + SAB_SECONDARY_SKEWED = 1 ! Skewed outgoing energy bins ! Elastic mode for S(a,b) elastic scattering integer, parameter :: & - & SAB_ELASTIC_DISCRETE = 3, & ! Sample from discrete cosines - & SAB_ELASTIC_EXACT = 4 ! Exact treatment for coherent elastic + SAB_ELASTIC_DISCRETE = 3, & ! Sample from discrete cosines + SAB_ELASTIC_EXACT = 4 ! Exact treatment for coherent elastic ! Reaction types integer, parameter :: & diff --git a/src/geometry.F90 b/src/geometry.F90 index 29ad60626c..0c5945c09b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -429,8 +429,8 @@ contains ! Couldn't find next cell anywhere! if ((.not. found) .and. (.not. plotting)) then message = "After particle crossed surface " // trim(int_to_str( & - surfaces(abs(p%surface)) % id)) // " it could not be located in " & - // "any cell and it did not leak." + surfaces(abs(p%surface)) % id)) // " it could not be located in & + &any cell and it did not leak." call fatal_error() end if @@ -1055,7 +1055,7 @@ contains y1 = surf % coeffs(5) z1 = surf % coeffs(6) if (x >= x0 .and. x < x1 .and. y >= y0 .and. y < y1 .and. & - & z >= z0 .and. z < z1) then + z >= z0 .and. z < z1) then s = SENSE_NEGATIVE else s = SENSE_POSITIVE @@ -1074,7 +1074,7 @@ contains I = surf % coeffs(9) J = surf % coeffs(10) func = A*x*x + B*y*y + C*z*z + D*x*y + E*y*z + F*x*z + G*x & - & + H*y + I*z + J + + H*y + I*z + J end select diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 5b3858755f..8fbf07353e 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -45,8 +45,8 @@ module geometry_header integer :: type ! Type of surface real(8), allocatable :: coeffs(:) ! Definition of surface integer, allocatable :: & - & neighbor_pos(:), & ! List of cells on positive side - & neighbor_neg(:) ! List of cells on negative side + neighbor_pos(:), & ! List of cells on positive side + neighbor_neg(:) ! List of cells on negative side integer :: bc ! Boundary condition end type Surface diff --git a/src/initialize.F90 b/src/initialize.F90 index 220c2a6809..bdcaa6b102 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -297,7 +297,7 @@ contains c % surfaces(j) = sign(i_array, id) else message = "Could not find surface " // trim(int_to_str(abs(id))) // & - & " specified on cell " // trim(int_to_str(c % id)) + " specified on cell " // trim(int_to_str(c % id)) call fatal_error() end if end if @@ -376,7 +376,7 @@ contains t % cell_bins(j) % scalar = dict_get_key(cell_dict, id) else message = "Could not find cell " // trim(int_to_str(id)) // & - & " specified on tally " // trim(int_to_str(t % id)) + " specified on tally " // trim(int_to_str(t % id)) call fatal_error() end if end do @@ -392,7 +392,7 @@ contains t % surface_bins(j) % scalar = dict_get_key(surface_dict, id) else message = "Could not find surface " // trim(int_to_str(id)) // & - & " specified on tally " // trim(int_to_str(t % id)) + " specified on tally " // trim(int_to_str(t % id)) call fatal_error() end if end do @@ -408,7 +408,7 @@ contains t % universe_bins(j) % scalar = dict_get_key(universe_dict, id) else message = "Could not find universe " // trim(int_to_str(id)) // & - & " specified on tally " // trim(int_to_str(t % id)) + " specified on tally " // trim(int_to_str(t % id)) call fatal_error() end if end do @@ -424,7 +424,7 @@ contains t % material_bins(j) % scalar = dict_get_key(material_dict, id) else message = "Could not find material " // trim(int_to_str(id)) // & - & " specified on tally " // trim(int_to_str(t % id)) + " specified on tally " // trim(int_to_str(t % id)) call fatal_error() end if end do @@ -440,7 +440,7 @@ contains t % cellborn_bins(j) % scalar = dict_get_key(cell_dict, id) else message = "Could not find material " // trim(int_to_str(id)) // & - & " specified on tally " // trim(int_to_str(t % id)) + " specified on tally " // trim(int_to_str(t % id)) call fatal_error() end if end do @@ -455,7 +455,7 @@ contains t % mesh = dict_get_key(mesh_dict, id) else message = "Could not find mesh " // trim(int_to_str(id)) // & - & " specified on tally " // trim(int_to_str(t % id)) + " specified on tally " // trim(int_to_str(t % id)) call fatal_error() end if end if @@ -545,9 +545,9 @@ contains ! Check to make sure either all atom percents or all weight percents are ! given if (.not. (all(mat%atom_percent > ZERO) .or. & - & all(mat%atom_percent < ZERO))) then + all(mat%atom_percent < ZERO))) then message = "Cannot mix atom and weight percents in material " // & - & int_to_str(mat % id) + int_to_str(mat % id) call fatal_error() end if @@ -605,7 +605,7 @@ contains end do sum_percent = ONE / sum_percent mat % density = -mat % density * N_AVOGADRO & - & / MASS_NEUTRON * sum_percent + / MASS_NEUTRON * sum_percent end if ! Calculate nuclide atom densities and deallocate atom_percent array diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e05edd524a..f96682724a 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -14,8 +14,8 @@ module input_xml implicit none - type(DictionaryII), pointer :: & ! used to count how many cells each - & cells_in_univ_dict => null() ! universe contains + type(DictionaryII), pointer :: & ! used to count how many cells each + cells_in_univ_dict => null() ! universe contains contains @@ -78,9 +78,8 @@ contains ! variable call get_environment_variable("CROSS_SECTIONS", env_variable) if (len_trim(env_variable) == 0) then - message = "No cross_sections.xml file was specified in " // & - "settings.xml or in the CROSS_SECTIONS environment " // & - "variable." + message = "No cross_sections.xml file was specified in settings.xml & + &or in the CROSS_SECTIONS environment variable." call fatal_error() else path_cross_sections = trim(env_variable) @@ -146,8 +145,8 @@ contains if (associated(entropy_box_)) then ! Check to make sure enough values were supplied if (size(entropy_box_) /= 6) then - message = "Need to supply lower-left and upper-right coordinates " & - // "for Shannon entropy box." + message = "Need to supply lower-left and upper-right coordinates & + &for Shannon entropy box." call fatal_error() end if @@ -157,8 +156,8 @@ contains ! Check on values provided if (.not. all(entropy_upper_right > entropy_lower_left)) then - message = "Upper-right coordinate must be greater than lower-left " & - // "coordinate for Shannon entropy box." + message = "Upper-right coordinate must be greater than lower-left & + &coordinate for Shannon entropy box." call fatal_error() end if @@ -336,7 +335,6 @@ contains if (n < coeffs_reqd) then message = "Not enough coefficients specified for surface: " // & trim(int_to_str(s % id)) - print *, n, coeffs_reqd call fatal_error() elseif (n > coeffs_reqd) then message = "Too many coefficients specified for surface: " // & @@ -562,8 +560,8 @@ contains trim(name) call fatal_error() elseif (nuc % ao /= ZERO .and. nuc % wo /= ZERO) then - message = "Cannot specify both atom and weight percents for a " & - // "nuclide: " // trim(name) + message = "Cannot specify both atom and weight percents for a & + &nuclide: " // trim(name) call fatal_error() end if @@ -693,16 +691,16 @@ contains ! Read mesh origin location if (m % n_dimension /= size(mesh_(i) % origin)) then - message = "Number of entries on must be the same as " // & - "the number of entries on ." + message = "Number of entries on must be the same as the & + &number of entries on ." call fatal_error() end if m % origin = mesh_(i) % origin ! Read mesh widths if (size(mesh_(i) % width) /= size(mesh_(i) % origin)) then - message = "Number of entries on must be the same as " // & - "the number of entries on ." + message = "Number of entries on must be the same as the & + &number of entries on ." call fatal_error() end if m % width = mesh_(i) % width @@ -845,8 +843,8 @@ contains case ('total') t % macro_bins(j) % scalar = MACRO_TOTAL if (t % n_bins(T_ENERGYOUT) > 0) then - message = "Cannot tally total reaction rate with an " & - // "outgoing energy filter." + message = "Cannot tally total reaction rate with an & + &outgoing energy filter." call fatal_error() end if case ('scatter') @@ -870,15 +868,15 @@ contains case ('absorption') t % macro_bins(j) % scalar = MACRO_ABSORPTION if (t % n_bins(T_ENERGYOUT) > 0) then - message = "Cannot tally absorption rate with an outgoing " & - // "energy filter." + message = "Cannot tally absorption rate with an outgoing & + &energy filter." call fatal_error() end if case ('fission') t % macro_bins(j) % scalar = MACRO_FISSION if (t % n_bins(T_ENERGYOUT) > 0) then - message = "Cannot tally fission rate with an outgoing " & - // "energy filter." + message = "Cannot tally fission rate with an outgoing & + &energy filter." call fatal_error() end if case ('nu-fission') @@ -890,9 +888,9 @@ contains ! Check to make sure that current is the only desired response ! for this tally if (n_words > 1) then - message = "Cannot tally other macro reactions in the same " & - // "tally as surface currents. Separate other macro " & - // "reactions into a distinct tally." + message = "Cannot tally other macro reactions in the same & + &tally as surface currents. Separate other macro & + &reactions into a distinct tally." call fatal_error() end if diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 3ec77bf6e0..c1405d743d 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -126,7 +126,7 @@ contains #ifdef MPI ! Collect number bank sites onto master process call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & - & MPI_COMM_WORLD, mpi_err) + MPI_COMM_WORLD, mpi_err) #else total_bank = n_bank #endif diff --git a/src/mpi_routines.F90 b/src/mpi_routines.F90 index 3729fc98f4..38029a36bd 100644 --- a/src/mpi_routines.F90 +++ b/src/mpi_routines.F90 @@ -82,7 +82,7 @@ contains bank_blocks = (/ 1, 3, 3, 1 /) bank_types = (/ MPI_INTEGER8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /) call MPI_TYPE_CREATE_STRUCT(4, bank_blocks, bank_disp, & - & bank_types, MPI_BANK, mpi_err) + bank_types, MPI_BANK, mpi_err) call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) #else @@ -115,9 +115,9 @@ contains integer(8) :: sites_needed ! # of sites to be sampled real(8) :: p_sample ! probability of sampling a site type(Bank), allocatable :: & - & temp_sites(:), & ! local array of extra sites on each node - & left_bank(:), & ! bank sites to send/recv to or from left node - & right_bank(:) ! bank sites to send/recv to or fram right node + temp_sites(:), & ! local array of extra sites on each node + left_bank(:), & ! bank sites to send/recv to or from left node + right_bank(:) ! bank sites to send/recv to or fram right node #ifdef MPI integer :: status(MPI_STATUS_SIZE) ! message status @@ -133,11 +133,11 @@ contains ! Determine starting index for fission bank and total sites in fission bank start = 0_8 call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, & - & MPI_COMM_WORLD, mpi_err) + MPI_COMM_WORLD, mpi_err) finish = start + n_bank total = finish call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, & - & MPI_COMM_WORLD, mpi_err) + MPI_COMM_WORLD, mpi_err) #else start = 0_8 @@ -200,11 +200,11 @@ contains #ifdef MPI start = 0_8 call MPI_EXSCAN(index_local, start, 1, MPI_INTEGER8, MPI_SUM, & - & MPI_COMM_WORLD, mpi_err) + MPI_COMM_WORLD, mpi_err) finish = start + index_local total = finish call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, & - & MPI_COMM_WORLD, mpi_err) + MPI_COMM_WORLD, mpi_err) #else start = 0_8 finish = index_local @@ -259,18 +259,18 @@ contains if (send_to_right > 0) then i = index_local - send_to_right + 1 call MPI_ISEND(temp_sites(i), send_to_right, MPI_BANK, rank+1, 0, & - & MPI_COMM_WORLD, request, mpi_err) + MPI_COMM_WORLD, request, mpi_err) else if (send_to_right < 0) then call MPI_IRECV(right_bank, -send_to_right, MPI_BANK, rank+1, 1, & - & MPI_COMM_WORLD, request_right, mpi_err) + MPI_COMM_WORLD, request_right, mpi_err) end if if (send_to_left < 0) then call MPI_IRECV(left_bank, -send_to_left, MPI_BANK, rank-1, 0, & - & MPI_COMM_WORLD, request_left, mpi_err) + MPI_COMM_WORLD, request_left, mpi_err) else if (send_to_left > 0) then call MPI_ISEND(temp_sites(1), send_to_left, MPI_BANK, rank-1, 1, & - & MPI_COMM_WORLD, request, mpi_err) + MPI_COMM_WORLD, request, mpi_err) end if #endif diff --git a/src/output.F90 b/src/output.F90 index 32e964123c..b09b97fb52 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -524,7 +524,7 @@ contains ! Write total atom density in atom/b-cm write(unit_,*) ' Atom Density = ' // trim(real_to_str(mat % density)) & - & // ' atom/b-cm' + // ' atom/b-cm' ! Write atom density for each nuclide in material write(unit_,*) ' Nuclides:' @@ -532,7 +532,7 @@ contains nuc => nuclides(mat % nuclide(i)) density = mat % atom_density(i) string = ' ' // trim(nuc % name) // ' = ' // & - & trim(real_to_str(density)) // ' atom/b-cm' + trim(real_to_str(density)) // ' atom/b-cm' write(unit_,*) trim(string) end do diff --git a/src/physics.F90 b/src/physics.F90 index 2a552c875b..449ee990b5 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -316,7 +316,7 @@ contains ! add to cumulative probability prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) & - & + f*(rxn%sigma(IE-rxn%IE+2))) + + f*(rxn%sigma(IE-rxn%IE+2))) end do else ! For nuclides with only total fission reaction, get a pointer to @@ -341,7 +341,7 @@ contains ! add to cumulative probability if (nuc % has_partial_fission) then prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) & - & + f*(rxn%sigma(IE-rxn%IE+2))) + + f*(rxn%sigma(IE-rxn%IE+2))) else prob = prob + micro_xs(index_nuclide) % fission end if @@ -437,7 +437,7 @@ contains ! add to cumulative probability prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) & - & + f*(rxn%sigma(IE-rxn%IE+2))) + + f*(rxn%sigma(IE-rxn%IE+2))) end do ! Perform collision physics for inelastics scattering @@ -1007,7 +1007,7 @@ contains ! determine outgoing energy in lab E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) & - & / ((A+ONE)*(A+ONE)) + / ((A+ONE)*(A+ONE)) ! determine outgoing angle in lab mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E) @@ -1083,7 +1083,7 @@ contains else i = binary_search(rxn % adist % energy, n, E) r = (E - rxn % adist % energy(i)) / & - & (rxn % adist % energy(i+1) - rxn % adist % energy(i)) + (rxn % adist % energy(i+1) - rxn % adist % energy(i)) end if ! Sample between the ith and (i+1)th bin @@ -1303,7 +1303,7 @@ contains lc = 2 + 2*NR i = binary_search(edist % data(lc+1), NE, E_in) r = (E_in - edist%data(lc+i)) / & - & (edist%data(lc+i+1) - edist%data(lc+i)) + (edist%data(lc+i+1) - edist%data(lc+i)) ! Sample outgoing energy bin r1 = prn() @@ -1379,7 +1379,7 @@ contains else i = binary_search(edist % data(lc+1), NE, E_in) r = (E_in - edist%data(lc+i)) / & - & (edist%data(lc+i+1) - edist%data(lc+i)) + (edist%data(lc+i+1) - edist%data(lc+i)) end if ! Sample between the ith and (i+1)th bin @@ -1452,7 +1452,7 @@ contains E_out = E_l_k + (r1 - c_k)/p_l_k else E_out = E_l_k + (sqrt(p_l_k*p_l_k + 2*frac*(r1 - c_k)) - & - & p_l_k)/frac + p_l_k)/frac end if else message = "Unknown interpolation type: " // trim(int_to_str(INTT)) @@ -1605,7 +1605,7 @@ contains else i = binary_search(edist % data(lc+1), NE, E_in) r = (E_in - edist%data(lc+i)) / & - & (edist%data(lc+i+1) - edist%data(lc+i)) + (edist%data(lc+i+1) - edist%data(lc+i)) end if ! Sample between the ith and (i+1)th bin @@ -1684,7 +1684,7 @@ contains E_out = E_l_k + (r1 - c_k)/p_l_k else E_out = E_l_k + (sqrt(p_l_k*p_l_k + 2*frac*(r1 - c_k)) - & - & p_l_k)/frac + p_l_k)/frac end if ! Determine Kalbach-Mann parameters @@ -1748,7 +1748,7 @@ contains else i = binary_search(edist % data(lc+1), NE, E_in) r = (E_in - edist%data(lc+i)) / & - & (edist%data(lc+i+1) - edist%data(lc+i)) + (edist%data(lc+i+1) - edist%data(lc+i)) end if ! Sample between the ith and (i+1)th bin @@ -1823,7 +1823,7 @@ contains E_out = E_l_k + (r1 - c_k)/p_l_k else E_out = E_l_k + (sqrt(p_l_k*p_l_k + 2*frac*(r1 - c_k)) - & - & p_l_k)/frac + p_l_k)/frac end if else message = "Unknown interpolation type: " // trim(int_to_str(INTT)) diff --git a/src/string.F90 b/src/string.F90 index 839abdea38..ae8991b849 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -348,7 +348,7 @@ contains ! Determine number of digits to right of decimal point index_decimal = index(string, '.') index_exponent = max(index(string, 'd'), index(string, 'D'), & - & index(string, 'e'), index(string, 'E')) + index(string, 'e'), index(string, 'E')) if (index_decimal > 0) then if (index_exponent > 0) then d = index_exponent - index_decimal - 1 From 2cf2e22bf8887c7b1b2602b13a5a8a92a4b10730 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Dec 2011 14:28:51 -0500 Subject: [PATCH 27/33] Removed fileio module. --- src/DEPENDENCIES | 5 --- src/OBJECTS | 1 - src/ace.F90 | 1 - src/fileio.F90 | 90 ------------------------------------------------ 4 files changed, 97 deletions(-) delete mode 100644 src/fileio.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 977875c1cb..bec099edaa 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -4,7 +4,6 @@ ace.o: datatypes.o ace.o: datatypes_header.o ace.o: endf.o ace.o: error.o -ace.o: fileio.o ace.o: fission.o ace.o: global.o ace.o: material_header.o @@ -39,10 +38,6 @@ energy_grid.o: datatypes_header.o energy_grid.o: global.o energy_grid.o: output.o -fileio.o: constants.o -fileio.o: global.o -fileio.o: string.o - fission.o: ace_header.o fission.o: constants.o fission.o: error.o diff --git a/src/OBJECTS b/src/OBJECTS index bbf3b2a257..caf0737868 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -10,7 +10,6 @@ endf.o \ endf_header.o \ energy_grid.o \ error.o \ -fileio.o \ fission.o \ geometry.o \ geometry_header.o \ diff --git a/src/ace.F90 b/src/ace.F90 index f4d5cf374e..3ab71869f4 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -8,7 +8,6 @@ module ace use datatypes_header, only: DictionaryCI, ListKeyValueCI use endf, only: reaction_name use error, only: fatal_error - use fileio, only: read_line, skip_lines use fission, only: nu_total use global use material_header, only: Material diff --git a/src/fileio.F90 b/src/fileio.F90 deleted file mode 100644 index 84aae2c83a..0000000000 --- a/src/fileio.F90 +++ /dev/null @@ -1,90 +0,0 @@ -module fileio - - use constants - use global, only: MAX_LINE_LEN, MAX_WORD_LEN, MAX_WORDS - use string, only: split_string_wl, lower_case - - implicit none - -contains - -!=============================================================================== -! READ_LINE reads a line from a file open on a unit -!=============================================================================== - - subroutine read_line(unit, line, ioError) - - integer, intent(in) :: unit ! unit to read from - character(*), intent(out) :: line ! line to return - integer, intent(out) :: ioError ! error status - - read(UNIT=unit, FMT='(A)', IOSTAT=ioError) line - - end subroutine read_line - -!=============================================================================== -! GET_NEXT_LINE reads the next line to the file connected on the specified unit -! including any continuation lines. If a line ends in an ampersand, the next -! line is read and its words are appended to the final array -!=============================================================================== - - subroutine get_next_line(unit, words, n, ioError) - - integer, intent(in) :: unit ! unit to read from - character(*), intent(out) :: words(MAX_WORDS) ! words read - integer, intent(out) :: n ! number of words - integer, intent(out) :: ioError ! error status - - character(MAX_LINE_LEN) :: line ! single line - character(MAX_WORD_LEN) :: local_words(MAX_WORDS) ! words on one line - integer :: index_word ! index of words - - index_word = 0 - do - ! read line from file - read(UNIT=unit, FMT='(A100)', IOSTAT=ioError) line - - ! if we're at the end of the file, return - if (ioError /= 0) return - - ! split a single line into words - call split_string_wl(line, local_words, n) - - ! if there are no words, we're done - if (n == 0) exit - - ! Check whether there is a continuation line - if (local_words(n) == '&') then - words(index_word+1:index_word+n-1) = local_words(1:n-1) - index_word = index_word + n - 1 - else - words(index_word+1:index_word+n) = local_words(1:n) - index_word = index_word + n - exit - end if - end do - - ! set total number of words - n = index_word - - end subroutine get_next_line - -!=============================================================================== -! SKIP_LINES skips 'n_lines' lines from a file open on a unit -!=============================================================================== - - subroutine skip_lines(unit, n_lines, ioError) - - integer, intent(in) :: unit ! unit to read from - integer, intent(in) :: n_lines ! number of lines to skip - integer, intent(out) :: ioError ! error status - - integer :: i ! index for number of lines - - do i = 1, n_lines - read(UNIT=unit, FMT=*, IOSTAT=ioError) - end do - - end subroutine skip_lines - -end module fileio From e25225d99334e3208924b25ca148e9b5f337e8ea Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Dec 2011 15:00:26 -0500 Subject: [PATCH 28/33] Moved MPI_BANK to global module. --- src/global.F90 | 1 + src/mpi_routines.F90 | 3 --- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/src/global.F90 b/src/global.F90 index e904f64a29..dc3d01092b 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -129,6 +129,7 @@ module global logical :: master ! master process? logical :: mpi_enabled ! is MPI in use and initialized? integer :: mpi_err ! MPI error code + integer :: MPI_BANK ! MPI datatype for fission bank ! ============================================================================ ! TIMING VARIABLES diff --git a/src/mpi_routines.F90 b/src/mpi_routines.F90 index 38029a36bd..8d6e07a076 100644 --- a/src/mpi_routines.F90 +++ b/src/mpi_routines.F90 @@ -15,9 +15,6 @@ module mpi_routines implicit none - integer :: MPI_BANK ! MPI datatype for fission bank - integer(8) :: bank_index ! Fission bank site unique identifier - contains !=============================================================================== From 06b0728aea16d4ab5478d1db450c7733bf9f7db6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Dec 2011 15:16:31 -0500 Subject: [PATCH 29/33] Moved mpi_routines subroutines into intercycle and initialize modules. --- src/DEPENDENCIES | 20 +-- src/OBJECTS | 1 - src/initialize.F90 | 81 ++++++++- src/intercycle.F90 | 322 +++++++++++++++++++++++++++++++++- src/main.F90 | 3 +- src/mpi_routines.F90 | 409 ------------------------------------------- src/tally.F90 | 2 +- 7 files changed, 409 insertions(+), 429 deletions(-) delete mode 100644 src/mpi_routines.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index bec099edaa..5016e4278b 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -68,6 +68,7 @@ global.o: tally_header.o global.o: timing.o initialize.o: ace.o +initialize.o: bank_header.o initialize.o: constants.o initialize.o: datatypes.o initialize.o: datatypes_header.o @@ -77,7 +78,6 @@ initialize.o: geometry.o initialize.o: geometry_header.o initialize.o: global.o initialize.o: input_xml.o -initialize.o: mpi_routines.o initialize.o: output.o initialize.o: random_lcg.o initialize.o: source.o @@ -100,9 +100,13 @@ input_xml.o: xml-fortran/templates/materials_t.o input_xml.o: xml-fortran/templates/settings_t.o input_xml.o: xml-fortran/templates/tallies_t.o -intercycle.o: global.o intercycle.o: error.o +intercycle.o: global.o intercycle.o: output.o +intercycle.o: particle_header.o +intercycle.o: random_lcg.o +intercycle.o: tally_header.o +intercycle.o: timing.o interpolation.o: constants.o interpolation.o: endf_header.o @@ -115,7 +119,6 @@ main.o: constants.o main.o: global.o main.o: initialize.o main.o: intercycle.o -main.o: mpi_routines.o main.o: output.o main.o: particle_header.o main.o: physics.o @@ -128,15 +131,6 @@ main.o: timing.o mesh.o: mesh_header.o -mpi_routines.o: constants.o -mpi_routines.o: error.o -mpi_routines.o: global.o -mpi_routines.o: output.o -mpi_routines.o: particle_header.o -mpi_routines.o: random_lcg.o -mpi_routines.o: tally_header.o -mpi_routines.o: timing.o - output.o: constants.o output.o: datatypes.o output.o: endf.o @@ -196,9 +190,9 @@ string.o: global.o tally.o: constants.o tally.o: error.o tally.o: global.o +tally.o: intercycle.o tally.o: mesh.o tally.o: mesh_header.o -tally.o: mpi_routines.o tally.o: output.o tally.o: search.o tally.o: string.o diff --git a/src/OBJECTS b/src/OBJECTS index caf0737868..e7fc0db9bd 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -22,7 +22,6 @@ main.o \ material_header.o \ mesh_header.o \ mesh.o \ -mpi_routines.o \ output.o \ particle_header.o \ physics.o \ diff --git a/src/initialize.F90 b/src/initialize.F90 index bdcaa6b102..248387a8b6 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -1,6 +1,7 @@ module initialize use ace, only: read_xs + use bank_header, only: Bank use constants use datatypes, only: dict_create, dict_add_key, dict_get_key, & dict_has_key, dict_keys @@ -12,7 +13,6 @@ module initialize use global use input_xml, only: read_input_xml, read_cross_sections_xml, & cells_in_univ_dict - use mpi_routines, only: setup_mpi use output, only: title, header, print_summary, print_geometry, & print_plot, create_summary_file use random_lcg, only: initialize_prng @@ -21,6 +21,10 @@ module initialize use tally, only: create_tally_map, TallyObject use timing, only: timer_start, timer_stop +#ifdef MPI + use mpi +#endif + implicit none type(DictionaryII), pointer :: build_dict => null() @@ -119,6 +123,81 @@ contains end subroutine initialize_run +!=============================================================================== +! SETUP_MPI initilizes the Message Passing Interface (MPI) and determines the +! number of processors the problem is being run with as well as the rank of each +! processor. +!=============================================================================== + + subroutine setup_mpi() + +#ifdef MPI + integer :: i + integer :: bank_blocks(4) ! Count for each datatype + integer :: bank_types(4) ! Datatypes + integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements + integer(MPI_ADDRESS_KIND) :: base + type(Bank) :: b + + mpi_enabled = .true. + + ! Initialize MPI + call MPI_INIT(mpi_err) + if (mpi_err /= MPI_SUCCESS) then + message = "Failed to initialize MPI." + call fatal_error() + end if + + ! Determine number of processors + call MPI_COMM_SIZE(MPI_COMM_WORLD, n_procs, mpi_err) + if (mpi_err /= MPI_SUCCESS) then + message = "Could not determine number of processors." + call fatal_error() + end if + + ! Determine rank of each processor + call MPI_COMM_RANK(MPI_COMM_WORLD, rank, mpi_err) + if (mpi_err /= MPI_SUCCESS) then + message = "Could not determine MPI rank." + call fatal_error() + end if + + ! Determine master + if (rank == 0) then + master = .true. + else + master = .false. + end if + + ! Determine displacements for MPI_BANK type + call MPI_GET_ADDRESS(b % id, bank_disp(1), mpi_err) + call MPI_GET_ADDRESS(b % xyz, bank_disp(2), mpi_err) + call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err) + call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err) + + ! Adjust displacements + base = bank_disp(1) + do i = 1, 4 + bank_disp(i) = bank_disp(i) - base + end do + + ! Define MPI_BANK for fission sites + bank_blocks = (/ 1, 3, 3, 1 /) + bank_types = (/ MPI_INTEGER8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /) + call MPI_TYPE_CREATE_STRUCT(4, bank_blocks, bank_disp, & + bank_types, MPI_BANK, mpi_err) + call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) + +#else + ! if no MPI, set processor to master + mpi_enabled = .false. + rank = 0 + n_procs = 1 + master = .true. +#endif + + end subroutine setup_mpi + !=============================================================================== ! READ_COMMAND_LINE reads all parameters from the command line !=============================================================================== diff --git a/src/intercycle.F90 b/src/intercycle.F90 index c1405d743d..50301927ce 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -2,9 +2,13 @@ module intercycle use ISO_FORTRAN_ENV + use error, only: fatal_error, warning use global - use error, only: warning - use output, only: write_message + use output, only: write_message + use particle_header, only: Particle, initialize_particle + use random_lcg, only: prn, set_particle_seed, prn_skip + use tally_header, only: TallyObject + use timing, only: timer_start, timer_stop #ifdef MPI use mpi @@ -12,6 +16,320 @@ module intercycle contains +!=============================================================================== +! SYNCHRONIZE_BANK samples source sites from the fission sites that were +! accumulated during the cycle. This routine is what allows this Monte Carlo to +! scale to large numbers of processors where other codes cannot. +!=============================================================================== + + subroutine synchronize_bank(i_cycle) + + integer, intent(in) :: i_cycle + + integer :: i, j, k ! loop indices + integer(8) :: start ! starting index in local fission bank + integer(8) :: finish ! ending index in local fission bank + integer(8) :: total ! total sites in global fission bank + integer(8) :: index_local ! index for source bank + integer :: send_to_left ! # of bank sites to send/recv to or from left + integer :: send_to_right ! # of bank sites to send/recv to or from right + integer(8) :: sites_needed ! # of sites to be sampled + real(8) :: p_sample ! probability of sampling a site + type(Bank), allocatable :: & + temp_sites(:), & ! local array of extra sites on each node + left_bank(:), & ! bank sites to send/recv to or from left node + right_bank(:) ! bank sites to send/recv to or fram right node + +#ifdef MPI + integer :: status(MPI_STATUS_SIZE) ! message status + integer :: request ! communication request for sending sites + integer :: request_left ! communication request for recv sites from left + integer :: request_right ! communication request for recv sites from right +#endif + + message = "Collecting number of fission sites..." + call write_message(8) + +#ifdef MPI + ! Determine starting index for fission bank and total sites in fission bank + start = 0_8 + call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, & + MPI_COMM_WORLD, mpi_err) + finish = start + n_bank + total = finish + call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, & + MPI_COMM_WORLD, mpi_err) + +#else + start = 0_8 + finish = n_bank + total = n_bank +#endif + + ! Check if there are no fission sites + if (total == 0) then + message = "No fission sites banked!" + call fatal_error() + end if + + ! Make sure all processors start at the same point for random sampling + call set_particle_seed(int(i_cycle,8)) + + ! Skip ahead however many random numbers are needed + call prn_skip(start) + + allocate(temp_sites(2*work)) + index_local = 0_8 ! Index for local source_bank + + if (total < n_particles) then + sites_needed = mod(n_particles,total) + else + sites_needed = n_particles + end if + p_sample = real(sites_needed,8)/real(total,8) + + message = "Sampling fission sites..." + call write_message(8) + + call timer_start(time_ic_sample) + + ! ========================================================================== + ! SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES + do i = 1, int(n_bank,4) + + ! If there are less than n_particles particles banked, automatically add + ! int(n_particles/total) sites to temp_sites. For example, if you need + ! 1000 and 300 were banked, this would add 3 source sites per banked site + ! and the remaining 100 would be randomly sampled. + if (total < n_particles) then + do j = 1,int(n_particles/total) + ! If index is within this node's range, add site to source + index_local = index_local + 1 + temp_sites(index_local) = fission_bank(i) + end do + end if + + ! Randomly sample sites needed + if (prn() < p_sample) then + index_local = index_local + 1 + temp_sites(index_local) = fission_bank(i) + end if + end do + + ! Now that we've sampled sites, check where the boundaries of data are for + ! the source bank +#ifdef MPI + start = 0_8 + call MPI_EXSCAN(index_local, start, 1, MPI_INTEGER8, MPI_SUM, & + MPI_COMM_WORLD, mpi_err) + finish = start + index_local + total = finish + call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, & + MPI_COMM_WORLD, mpi_err) +#else + start = 0_8 + finish = index_local + total = index_local +#endif + + ! Determine how many sites to send to adjacent nodes + send_to_left = int(bank_first - 1_8 - start, 4) + send_to_right = int(finish - bank_last, 4) + + ! Check to make sure number of sites is not more than size of bank + if (abs(send_to_left) > work .or. abs(send_to_right) > work) then + message = "Tried sending sites to neighboring process greater than " & + // "the size of the source bank." + call fatal_error() + end if + + if (rank == n_procs - 1) then + if (total > n_particles) then + ! If we have extra sites sampled, we will simply discard the extra + ! ones on the last processor + if (rank == n_procs - 1) then + index_local = index_local - send_to_right + end if + + elseif (total < n_particles) then + ! If we have too few sites, grab sites from the very end of the + ! fission bank + sites_needed = n_particles - total + do i = 1, int(sites_needed,4) + index_local = index_local + 1 + temp_sites(index_local) = fission_bank(n_bank - sites_needed + i) + end do + end if + + ! the last processor should not be sending sites to right + send_to_right = 0 + end if + + call timer_stop(time_ic_sample) + call timer_start(time_ic_sendrecv) + +#ifdef MPI + message = "Sending fission sites..." + call write_message(8) + + ! ========================================================================== + ! SEND BANK SITES TO NEIGHBORS + allocate(left_bank(abs(send_to_left))) + allocate(right_bank(abs(send_to_right))) + + if (send_to_right > 0) then + i = index_local - send_to_right + 1 + call MPI_ISEND(temp_sites(i), send_to_right, MPI_BANK, rank+1, 0, & + MPI_COMM_WORLD, request, mpi_err) + else if (send_to_right < 0) then + call MPI_IRECV(right_bank, -send_to_right, MPI_BANK, rank+1, 1, & + MPI_COMM_WORLD, request_right, mpi_err) + end if + + if (send_to_left < 0) then + call MPI_IRECV(left_bank, -send_to_left, MPI_BANK, rank-1, 0, & + MPI_COMM_WORLD, request_left, mpi_err) + else if (send_to_left > 0) then + call MPI_ISEND(temp_sites(1), send_to_left, MPI_BANK, rank-1, 1, & + MPI_COMM_WORLD, request, mpi_err) + end if +#endif + + call timer_stop(time_ic_sendrecv) + call timer_start(time_ic_rebuild) + + message = "Constructing source bank..." + call write_message(8) + + ! ========================================================================== + ! RECONSTRUCT SOURCE BANK + if (send_to_left < 0 .and. send_to_right >= 0) then + i = -send_to_left ! size of first block + j = int(index_local,4) - send_to_right ! size of second block + call copy_from_bank(temp_sites, i+1, j) +#ifdef MPI + call MPI_WAIT(request_left, status, mpi_err) +#endif + call copy_from_bank(left_bank, 1, i) + else if (send_to_left >= 0 .and. send_to_right < 0) then + i = int(index_local,4) - send_to_left ! size of first block + j = -send_to_right ! size of second block + call copy_from_bank(temp_sites(1+send_to_left), 1, i) +#ifdef MPI + call MPI_WAIT(request_right, status, mpi_err) +#endif + call copy_from_bank(right_bank, i+1, j) + else if (send_to_left >= 0 .and. send_to_right >= 0) then + i = int(index_local,4) - send_to_left - send_to_right + call copy_from_bank(temp_sites(1+send_to_left), 1, i) + else if (send_to_left < 0 .and. send_to_right < 0) then + i = -send_to_left + j = int(index_local,4) + k = -send_to_right + call copy_from_bank(temp_sites, i+1, j) +#ifdef MPI + call MPI_WAIT(request_left, status, mpi_err) +#endif + call copy_from_bank(left_bank, 1, i) +#ifdef MPI + call MPI_WAIT(request_right, status, mpi_err) +#endif + call copy_from_bank(right_bank, i+j+1, k) + end if + + ! Reset source index + source_index = 0_8 + + call timer_stop(time_ic_rebuild) + +#ifdef MPI + deallocate(left_bank) + deallocate(right_bank) +#endif + deallocate(temp_sites) + + end subroutine synchronize_bank + +!=============================================================================== +! COPY_FROM_BANK +!=============================================================================== + + subroutine copy_from_bank(temp_bank, i_start, n_sites) + + integer, intent(in) :: n_sites ! # of bank sites to copy + type(Bank), intent(in) :: temp_bank(n_sites) + integer, intent(in) :: i_start ! starting index in source_bank + + integer :: i ! index in temp_bank + integer :: i_source ! index in source_bank + type(Particle), pointer :: p + + do i = 1, n_sites + i_source = i_start + i - 1 + p => source_bank(i_source) + + ! set defaults + call initialize_particle(p) + + p % coord % xyz = temp_bank(i) % xyz + p % coord % uvw = temp_bank(i) % uvw + p % last_xyz = temp_bank(i) % xyz + p % E = temp_bank(i) % E + p % last_E = temp_bank(i) % E + + end do + + end subroutine copy_from_bank + +!=============================================================================== +! REDUCE_TALLIES collects all the results from tallies onto one processor +!=============================================================================== + +#ifdef MPI + subroutine reduce_tallies() + + integer :: i + integer :: n + integer :: m + integer :: n_bins + real(8), allocatable :: tally_temp(:,:) + type(TallyObject), pointer :: t + + do i = 1, n_tallies + t => tallies(i) + + n = t % n_total_bins + m = t % n_macro_bins + n_bins = n*m + + allocate(tally_temp(n,m)) + + tally_temp = t % scores(:,:) % val_history + + if (master) then + ! The MPI_IN_PLACE specifier allows the master to copy values into a + ! receive buffer without having a temporary variable + call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, MPI_SUM, & + 0, MPI_COMM_WORLD, mpi_err) + + ! Transfer values to val_history on master + t % scores(:,:) % val_history = tally_temp + else + ! Receive buffer not significant at other processors + call MPI_REDUCE(tally_temp, tally_temp, n_bins, MPI_REAL8, MPI_SUM, & + 0, MPI_COMM_WORLD, mpi_err) + + ! Reset val_history on other processors + t % scores(:,:) % val_history = 0 + end if + + deallocate(tally_temp) + + end do + + end subroutine reduce_tallies +#endif + !=============================================================================== ! SHANNON_ENTROPY calculates the Shannon entropy of the fission source ! distribution to assess source convergence diff --git a/src/main.F90 b/src/main.F90 index ca0094b140..d71993c923 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -3,8 +3,7 @@ program main use constants use global use initialize, only: initialize_run - use intercycle, only: shannon_entropy, calculate_keff - use mpi_routines, only: synchronize_bank + use intercycle, only: shannon_entropy, calculate_keff, synchronize_bank use output, only: write_message, header, print_runtime use particle_header, only: Particle use plot, only: run_plot diff --git a/src/mpi_routines.F90 b/src/mpi_routines.F90 deleted file mode 100644 index 8d6e07a076..0000000000 --- a/src/mpi_routines.F90 +++ /dev/null @@ -1,409 +0,0 @@ -module mpi_routines - - use constants, only: MAX_LINE_LEN - use error, only: fatal_error - use global - use output, only: write_message - use particle_header, only: Particle, initialize_particle - use random_lcg, only: prn, set_particle_seed, prn_skip - use tally_header, only: TallyObject - use timing, only: timer_start, timer_stop - -#ifdef MPI - use mpi -#endif - - implicit none - -contains - -!=============================================================================== -! SETUP_MPI initilizes the Message Passing Interface (MPI) and determines the -! number of processors the problem is being run with as well as the rank of each -! processor. -!=============================================================================== - - subroutine setup_mpi() - -#ifdef MPI - integer :: i - integer :: bank_blocks(4) ! Count for each datatype - integer :: bank_types(4) ! Datatypes - integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements - integer(MPI_ADDRESS_KIND) :: base - type(Bank) :: b - - mpi_enabled = .true. - - ! Initialize MPI - call MPI_INIT(mpi_err) - if (mpi_err /= MPI_SUCCESS) then - message = "Failed to initialize MPI." - call fatal_error() - end if - - ! Determine number of processors - call MPI_COMM_SIZE(MPI_COMM_WORLD, n_procs, mpi_err) - if (mpi_err /= MPI_SUCCESS) then - message = "Could not determine number of processors." - call fatal_error() - end if - - ! Determine rank of each processor - call MPI_COMM_RANK(MPI_COMM_WORLD, rank, mpi_err) - if (mpi_err /= MPI_SUCCESS) then - message = "Could not determine MPI rank." - call fatal_error() - end if - - ! Determine master - if (rank == 0) then - master = .true. - else - master = .false. - end if - - ! Determine displacements for MPI_BANK type - call MPI_GET_ADDRESS(b % id, bank_disp(1), mpi_err) - call MPI_GET_ADDRESS(b % xyz, bank_disp(2), mpi_err) - call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err) - call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err) - - ! Adjust displacements - base = bank_disp(1) - do i = 1, 4 - bank_disp(i) = bank_disp(i) - base - end do - - ! Define MPI_BANK for fission sites - bank_blocks = (/ 1, 3, 3, 1 /) - bank_types = (/ MPI_INTEGER8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /) - call MPI_TYPE_CREATE_STRUCT(4, bank_blocks, bank_disp, & - bank_types, MPI_BANK, mpi_err) - call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) - -#else - ! if no MPI, set processor to master - mpi_enabled = .false. - rank = 0 - n_procs = 1 - master = .true. -#endif - - end subroutine setup_mpi - -!=============================================================================== -! SYNCHRONIZE_BANK samples source sites from the fission sites that were -! accumulated during the cycle. This routine is what allows this Monte Carlo to -! scale to large numbers of processors where other codes cannot. -!=============================================================================== - - subroutine synchronize_bank(i_cycle) - - integer, intent(in) :: i_cycle - - integer :: i, j, k ! loop indices - integer(8) :: start ! starting index in local fission bank - integer(8) :: finish ! ending index in local fission bank - integer(8) :: total ! total sites in global fission bank - integer(8) :: index_local ! index for source bank - integer :: send_to_left ! # of bank sites to send/recv to or from left - integer :: send_to_right ! # of bank sites to send/recv to or from right - integer(8) :: sites_needed ! # of sites to be sampled - real(8) :: p_sample ! probability of sampling a site - type(Bank), allocatable :: & - temp_sites(:), & ! local array of extra sites on each node - left_bank(:), & ! bank sites to send/recv to or from left node - right_bank(:) ! bank sites to send/recv to or fram right node - -#ifdef MPI - integer :: status(MPI_STATUS_SIZE) ! message status - integer :: request ! communication request for sending sites - integer :: request_left ! communication request for recv sites from left - integer :: request_right ! communication request for recv sites from right -#endif - - message = "Collecting number of fission sites..." - call write_message(8) - -#ifdef MPI - ! Determine starting index for fission bank and total sites in fission bank - start = 0_8 - call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, & - MPI_COMM_WORLD, mpi_err) - finish = start + n_bank - total = finish - call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, & - MPI_COMM_WORLD, mpi_err) - -#else - start = 0_8 - finish = n_bank - total = n_bank -#endif - - ! Check if there are no fission sites - if (total == 0) then - message = "No fission sites banked!" - call fatal_error() - end if - - ! Make sure all processors start at the same point for random sampling - call set_particle_seed(int(i_cycle,8)) - - ! Skip ahead however many random numbers are needed - call prn_skip(start) - - allocate(temp_sites(2*work)) - index_local = 0_8 ! Index for local source_bank - - if (total < n_particles) then - sites_needed = mod(n_particles,total) - else - sites_needed = n_particles - end if - p_sample = real(sites_needed,8)/real(total,8) - - message = "Sampling fission sites..." - call write_message(8) - - call timer_start(time_ic_sample) - - ! ========================================================================== - ! SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES - do i = 1, int(n_bank,4) - - ! If there are less than n_particles particles banked, automatically add - ! int(n_particles/total) sites to temp_sites. For example, if you need - ! 1000 and 300 were banked, this would add 3 source sites per banked site - ! and the remaining 100 would be randomly sampled. - if (total < n_particles) then - do j = 1,int(n_particles/total) - ! If index is within this node's range, add site to source - index_local = index_local + 1 - temp_sites(index_local) = fission_bank(i) - end do - end if - - ! Randomly sample sites needed - if (prn() < p_sample) then - index_local = index_local + 1 - temp_sites(index_local) = fission_bank(i) - end if - end do - - ! Now that we've sampled sites, check where the boundaries of data are for - ! the source bank -#ifdef MPI - start = 0_8 - call MPI_EXSCAN(index_local, start, 1, MPI_INTEGER8, MPI_SUM, & - MPI_COMM_WORLD, mpi_err) - finish = start + index_local - total = finish - call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, & - MPI_COMM_WORLD, mpi_err) -#else - start = 0_8 - finish = index_local - total = index_local -#endif - - ! Determine how many sites to send to adjacent nodes - send_to_left = int(bank_first - 1_8 - start, 4) - send_to_right = int(finish - bank_last, 4) - - ! Check to make sure number of sites is not more than size of bank - if (abs(send_to_left) > work .or. abs(send_to_right) > work) then - message = "Tried sending sites to neighboring process greater than " & - // "the size of the source bank." - call fatal_error() - end if - - if (rank == n_procs - 1) then - if (total > n_particles) then - ! If we have extra sites sampled, we will simply discard the extra - ! ones on the last processor - if (rank == n_procs - 1) then - index_local = index_local - send_to_right - end if - - elseif (total < n_particles) then - ! If we have too few sites, grab sites from the very end of the - ! fission bank - sites_needed = n_particles - total - do i = 1, int(sites_needed,4) - index_local = index_local + 1 - temp_sites(index_local) = fission_bank(n_bank - sites_needed + i) - end do - end if - - ! the last processor should not be sending sites to right - send_to_right = 0 - end if - - call timer_stop(time_ic_sample) - call timer_start(time_ic_sendrecv) - -#ifdef MPI - message = "Sending fission sites..." - call write_message(8) - - ! ========================================================================== - ! SEND BANK SITES TO NEIGHBORS - allocate(left_bank(abs(send_to_left))) - allocate(right_bank(abs(send_to_right))) - - if (send_to_right > 0) then - i = index_local - send_to_right + 1 - call MPI_ISEND(temp_sites(i), send_to_right, MPI_BANK, rank+1, 0, & - MPI_COMM_WORLD, request, mpi_err) - else if (send_to_right < 0) then - call MPI_IRECV(right_bank, -send_to_right, MPI_BANK, rank+1, 1, & - MPI_COMM_WORLD, request_right, mpi_err) - end if - - if (send_to_left < 0) then - call MPI_IRECV(left_bank, -send_to_left, MPI_BANK, rank-1, 0, & - MPI_COMM_WORLD, request_left, mpi_err) - else if (send_to_left > 0) then - call MPI_ISEND(temp_sites(1), send_to_left, MPI_BANK, rank-1, 1, & - MPI_COMM_WORLD, request, mpi_err) - end if -#endif - - call timer_stop(time_ic_sendrecv) - call timer_start(time_ic_rebuild) - - message = "Constructing source bank..." - call write_message(8) - - ! ========================================================================== - ! RECONSTRUCT SOURCE BANK - if (send_to_left < 0 .and. send_to_right >= 0) then - i = -send_to_left ! size of first block - j = int(index_local,4) - send_to_right ! size of second block - call copy_from_bank(temp_sites, i+1, j) -#ifdef MPI - call MPI_WAIT(request_left, status, mpi_err) -#endif - call copy_from_bank(left_bank, 1, i) - else if (send_to_left >= 0 .and. send_to_right < 0) then - i = int(index_local,4) - send_to_left ! size of first block - j = -send_to_right ! size of second block - call copy_from_bank(temp_sites(1+send_to_left), 1, i) -#ifdef MPI - call MPI_WAIT(request_right, status, mpi_err) -#endif - call copy_from_bank(right_bank, i+1, j) - else if (send_to_left >= 0 .and. send_to_right >= 0) then - i = int(index_local,4) - send_to_left - send_to_right - call copy_from_bank(temp_sites(1+send_to_left), 1, i) - else if (send_to_left < 0 .and. send_to_right < 0) then - i = -send_to_left - j = int(index_local,4) - k = -send_to_right - call copy_from_bank(temp_sites, i+1, j) -#ifdef MPI - call MPI_WAIT(request_left, status, mpi_err) -#endif - call copy_from_bank(left_bank, 1, i) -#ifdef MPI - call MPI_WAIT(request_right, status, mpi_err) -#endif - call copy_from_bank(right_bank, i+j+1, k) - end if - - ! Reset source index - source_index = 0_8 - - call timer_stop(time_ic_rebuild) - -#ifdef MPI - deallocate(left_bank) - deallocate(right_bank) -#endif - deallocate(temp_sites) - - end subroutine synchronize_bank - -!=============================================================================== -! COPY_FROM_BANK -!=============================================================================== - - subroutine copy_from_bank(temp_bank, i_start, n_sites) - - integer, intent(in) :: n_sites ! # of bank sites to copy - type(Bank), intent(in) :: temp_bank(n_sites) - integer, intent(in) :: i_start ! starting index in source_bank - - integer :: i ! index in temp_bank - integer :: i_source ! index in source_bank - type(Particle), pointer :: p - - do i = 1, n_sites - i_source = i_start + i - 1 - p => source_bank(i_source) - - ! set defaults - call initialize_particle(p) - - p % coord % xyz = temp_bank(i) % xyz - p % coord % uvw = temp_bank(i) % uvw - p % last_xyz = temp_bank(i) % xyz - p % E = temp_bank(i) % E - p % last_E = temp_bank(i) % E - - end do - - end subroutine copy_from_bank - -!=============================================================================== -! REDUCE_TALLIES collects all the results from tallies onto one processor -!=============================================================================== - -#ifdef MPI - subroutine reduce_tallies() - - integer :: i - integer :: n - integer :: m - integer :: n_bins - real(8), allocatable :: tally_temp(:,:) - type(TallyObject), pointer :: t - - do i = 1, n_tallies - t => tallies(i) - - n = t % n_total_bins - m = t % n_macro_bins - n_bins = n*m - - allocate(tally_temp(n,m)) - - tally_temp = t % scores(:,:) % val_history - - if (master) then - ! The MPI_IN_PLACE specifier allows the master to copy values into a - ! receive buffer without having a temporary variable - call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) - - ! Transfer values to val_history on master - t % scores(:,:) % val_history = tally_temp - else - ! Receive buffer not significant at other processors - call MPI_REDUCE(tally_temp, tally_temp, n_bins, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) - - ! Reset val_history on other processors - t % scores(:,:) % val_history = 0 - end if - - deallocate(tally_temp) - - end do - - end subroutine reduce_tallies -#endif - -end module mpi_routines diff --git a/src/tally.F90 b/src/tally.F90 index 8b039c8af6..23d0e9bc9e 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -11,8 +11,8 @@ module tally use tally_header, only: TallyScore, TallyMapItem, TallyMapElement #ifdef MPI + use intercycle, only: reduce_tallies use mpi - use mpi_routines, only: reduce_tallies #endif implicit none From bcebb0ccd1f6af7661169a4ea337586b4f02ef78 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Dec 2011 15:48:16 -0500 Subject: [PATCH 30/33] Enhanced fatal_error subroutine. --- src/error.F90 | 53 +++++++++++++++++++++++++++++++++----------------- src/global.F90 | 5 ----- src/main.F90 | 5 +++++ 3 files changed, 40 insertions(+), 23 deletions(-) diff --git a/src/error.F90 b/src/error.F90 index fab7c8dc84..2eac3f2372 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -2,7 +2,11 @@ module error use ISO_FORTRAN_ENV - use global, only: master, free_memory, message + use global, only: master, free_memory, message, mpi_err + +#ifdef MPI + use mpi +#endif implicit none @@ -25,14 +29,14 @@ contains ! Only allow master to print to screen if (.not. master) return - write(ou, fmt='(1X,A9)', advance='no') 'WARNING: ' + write(OUTPUT_UNIT, fmt='(1X,A9)', advance='no') 'WARNING: ' n_lines = (len_trim(message)-1)/70 + 1 do i = 1, n_lines if (i == 1) then - write(ou, fmt='(A70)') message(70*(i-1)+1:70*i) + write(OUTPUT_UNIT, fmt='(A70)') message(70*(i-1)+1:70*i) else - write(ou, fmt='(10X,A70)') message(70*(i-1)+1:70*i) + write(OUTPUT_UNIT, fmt='(10X,A70)') message(70*(i-1)+1:70*i) end if end do @@ -44,29 +48,42 @@ contains ! the program is aborted. !=============================================================================== - subroutine fatal_error() + subroutine fatal_error(error_code) + integer, optional :: error_code ! error code + + integer :: code ! error code integer :: n_lines ! number of lines integer :: i ! loop index over lines - ! Only allow master to print to screen - if (master) then - write(eu, fmt='(1X,A7)', advance='no') 'ERROR: ' - - n_lines = (len_trim(message)-1)/72 + 1 - do i = 1, n_lines - if (i == 1) then - write(eu, fmt='(A72)') message(72*(i-1)+1:72*i) - else - write(eu, fmt='(7X,A72)') message(72*(i-1)+1:72*i) - end if - end do - write(eu,*) + ! set default error code + if (present(error_code)) then + code = error_code + else + code = -1 end if + ! Only allow master to print to screen + write(ERROR_UNIT, fmt='(1X,A7)', advance='no') 'ERROR: ' + + n_lines = (len_trim(message)-1)/72 + 1 + do i = 1, n_lines + if (i == 1) then + write(ERROR_UNIT, fmt='(A72)') message(72*(i-1)+1:72*i) + else + write(ERROR_UNIT, fmt='(7X,A72)') message(72*(i-1)+1:72*i) + end if + end do + write(ERROR_UNIT,*) + ! Release memory from all allocatable arrays call free_memory() +#ifdef MPI + ! Abort MPI + call MPI_ABORT(MPI_COMM_WORLD, code, mpi_err) +#endif + ! Abort program stop diff --git a/src/global.F90 b/src/global.F90 index dc3d01092b..347eace57b 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -211,11 +211,6 @@ contains if (allocated(fission_bank)) deallocate(fission_bank) if (allocated(source_bank)) deallocate(source_bank) -#ifdef MPI - ! If MPI is in use and enabled, terminate it - call MPI_FINALIZE(mpi_err) -#endif - end subroutine free_memory end module global diff --git a/src/main.F90 b/src/main.F90 index d71993c923..9586f8bc40 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -44,6 +44,11 @@ program main ! deallocate arrays call free_memory() + +#ifdef MPI + ! If MPI is in use and enabled, terminate it + call MPI_FINALIZE(mpi_err) +#endif contains From fd9cc4019380d2836be4baa329d5e7252cc1a389 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Dec 2011 23:28:13 -0500 Subject: [PATCH 31/33] Changed module interface int_to_str to to_str. --- src/ace.F90 | 4 +- src/endf.F90 | 16 ++-- src/geometry.F90 | 34 ++++----- src/initialize.F90 | 48 ++++++------ src/input_xml.F90 | 26 +++---- src/interpolation.F90 | 6 +- src/main.F90 | 4 +- src/output.F90 | 168 +++++++++++++++++++++--------------------- src/physics.F90 | 24 +++--- src/source.F90 | 8 +- src/string.F90 | 4 +- src/tally.F90 | 84 ++++++++++----------- 12 files changed, 213 insertions(+), 213 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index 3ab71869f4..32894116c3 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -13,7 +13,7 @@ module ace use material_header, only: Material use output, only: write_message, print_nuclide, header use string, only: split_string, str_to_int, str_to_real, & - lower_case, int_to_str + lower_case, to_str implicit none @@ -203,7 +203,7 @@ contains ! Check to make sure S(a,b) table matched a nuclide if (mat % sab_nuclide == 0) then message = "S(a,b) table " // trim(mat % sab_name) // " did not match & - &any nuclide on material " // trim(int_to_str(mat % id)) + &any nuclide on material " // trim(to_str(mat % id)) call fatal_error() end if end if diff --git a/src/endf.F90 b/src/endf.F90 index 5f35da5b07..92103b5c7e 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -1,7 +1,7 @@ module endf use constants - use string, only: int_to_str + use string, only: to_str contains @@ -72,7 +72,7 @@ contains case (N_NPA) string = '(n,npa)' case (N_N1 : N_N40) - string = '(n,n' // trim(int_to_str(MT-50)) // ')' + string = '(n,n' // trim(to_str(MT-50)) // ')' case (N_NC) string = '(n,nc)' case (N_GAMMA) @@ -122,27 +122,27 @@ contains case (444) string = '(damage)' case (600 : 648) - string = '(n,p' // trim(int_to_str(MT-600)) // ')' + string = '(n,p' // trim(to_str(MT-600)) // ')' case (649) string = '(n,pc)' case (650 : 698) - string = '(n,d' // trim(int_to_str(MT-650)) // ')' + string = '(n,d' // trim(to_str(MT-650)) // ')' case (699) string = '(n,dc)' case (700 : 748) - string = '(n,t' // trim(int_to_str(MT-700)) // ')' + string = '(n,t' // trim(to_str(MT-700)) // ')' case (749) string = '(n,tc)' case (750 : 798) - string = '(n,3He' // trim(int_to_str(MT-750)) // ')' + string = '(n,3He' // trim(to_str(MT-750)) // ')' case (799) string = '(n,3Hec)' case (800 : 848) - string = '(n,a' // trim(int_to_str(MT-800)) // ')' + string = '(n,a' // trim(to_str(MT-800)) // ')' case (849) string = '(n,tc)' case default - string = 'MT=' // trim(int_to_str(MT)) + string = 'MT=' // trim(to_str(MT)) end select end function reaction_name diff --git a/src/geometry.F90 b/src/geometry.F90 index 0c5945c09b..983cf5d3c1 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -7,7 +7,7 @@ module geometry use global use output, only: write_message use particle_header, only: Particle, LocalCoord, deallocate_coord - use string, only: int_to_str + use string, only: to_str use tally, only: score_surface_current implicit none @@ -249,7 +249,7 @@ contains surf => surfaces(abs(p % surface)) if (verbosity >= 10 .or. trace) then - message = " Crossing surface " // trim(int_to_str(surf % id)) + message = " Crossing surface " // trim(to_str(surf % id)) call write_message() end if @@ -274,7 +274,7 @@ contains ! Display message if (verbosity >= 10 .or. trace) then - message = " Leaked out of surface " // trim(int_to_str(surf % id)) + message = " Leaked out of surface " // trim(to_str(surf % id)) call write_message() end if return @@ -383,7 +383,7 @@ contains p % coord0 % uvw = (/ u, v, w /) case default message = "Reflection not supported for surface " // & - trim(int_to_str(surf % id)) + trim(to_str(surf % id)) call fatal_error() end select @@ -396,7 +396,7 @@ contains ! Diagnostic message if (verbosity >= 10 .or. trace) then - message = " Reflected from surface " // trim(int_to_str(surf%id)) + message = " Reflected from surface " // trim(to_str(surf%id)) call write_message() end if return @@ -428,7 +428,7 @@ contains ! Couldn't find next cell anywhere! if ((.not. found) .and. (.not. plotting)) then - message = "After particle crossed surface " // trim(int_to_str( & + message = "After particle crossed surface " // trim(to_str( & surfaces(abs(p%surface)) % id)) // " it could not be located in & &any cell and it did not leak." call fatal_error() @@ -455,9 +455,9 @@ contains lat => lattices(p % coord % lattice) if (verbosity >= 10 .or. trace) then - message = " Crossing lattice " // trim(int_to_str(lat % id)) // & - ". Current position (" // trim(int_to_str(p % coord % lattice_x)) & - // "," // trim(int_to_str(p % coord % lattice_y)) // ")" + message = " Crossing lattice " // trim(to_str(lat % id)) // & + ". Current position (" // trim(to_str(p % coord % lattice_x)) & + // "," // trim(to_str(p % coord % lattice_y)) // ")" call write_message() end if @@ -495,14 +495,14 @@ contains i_x = p % coord % lattice_x i_y = p % coord % lattice_y if (i_x < 1 .or. i_x > lat % n_x) then - message = "Reached edge of lattice " // trim(int_to_str(lat % id)) // & - " at position (" // trim(int_to_str(i_x)) // "," // & - trim(int_to_str(i_y)) // ")." + message = "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 = "Reached edge of lattice " // trim(int_to_str(lat % id)) // & - " at position (" // trim(int_to_str(i_x)) // "," // & - trim(int_to_str(i_y)) // ")." + message = "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 @@ -512,8 +512,8 @@ contains ! Find cell in next lattice element call find_cell(p, found) if (.not. found) then - message = "Could not locate particle in universe: " // & - int_to_str(universes(p % coord % universe) % id) + message = "Could not locate particle " // trim(to_str(p % id)) // & + " in universe " // to_str(universes(p % coord % universe) % id) call fatal_error() end if diff --git a/src/initialize.F90 b/src/initialize.F90 index 248387a8b6..1bcae2ef36 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -17,7 +17,7 @@ module initialize print_plot, create_summary_file use random_lcg, only: initialize_prng use source, only: initialize_source - use string, only: int_to_str, starts_with, ends_with, lower_case + use string, only: to_str, starts_with, ends_with, lower_case use tally, only: create_tally_map, TallyObject use timing, only: timer_start, timer_stop @@ -375,8 +375,8 @@ contains i_array = dict_get_key(surface_dict, abs(id)) c % surfaces(j) = sign(i_array, id) else - message = "Could not find surface " // trim(int_to_str(abs(id))) // & - " specified on cell " // trim(int_to_str(c % id)) + message = "Could not find surface " // trim(to_str(abs(id))) // & + " specified on cell " // trim(to_str(c % id)) call fatal_error() end if end if @@ -389,8 +389,8 @@ contains if (dict_has_key(universe_dict, id)) then c % universe = dict_get_key(universe_dict, id) else - message = "Could not find universe " // trim(int_to_str(id)) // & - " specified on cell " // trim(int_to_str(c % id)) + message = "Could not find universe " // trim(to_str(id)) // & + " specified on cell " // trim(to_str(c % id)) call fatal_error() end if @@ -403,8 +403,8 @@ contains c % type = CELL_NORMAL c % material = dict_get_key(material_dict, id) else - message = "Could not find material " // trim(int_to_str(id)) // & - " specified on cell " // trim(int_to_str(c % id)) + message = "Could not find material " // trim(to_str(id)) // & + " specified on cell " // trim(to_str(c % id)) call fatal_error() end if else @@ -416,8 +416,8 @@ contains c % type = CELL_LATTICE c % fill = dict_get_key(lattice_dict, id) else - message = "Specified fill " // trim(int_to_str(id)) // " on cell " // & - trim(int_to_str(c % id)) // " is neither a universe nor a lattice." + message = "Specified fill " // trim(to_str(id)) // " on cell " // & + trim(to_str(c % id)) // " is neither a universe nor a lattice." call fatal_error() end if end if @@ -434,8 +434,8 @@ contains if (dict_has_key(universe_dict, id)) then l % element(j,k) = dict_get_key(universe_dict, id) else - message = "Invalid universe number " // trim(int_to_str(id)) & - // " specified on lattice " // trim(int_to_str(l % id)) + message = "Invalid universe number " // trim(to_str(id)) & + // " specified on lattice " // trim(to_str(l % id)) call fatal_error() end if end do @@ -454,8 +454,8 @@ contains if (dict_has_key(cell_dict, id)) then t % cell_bins(j) % scalar = dict_get_key(cell_dict, id) else - message = "Could not find cell " // trim(int_to_str(id)) // & - " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find cell " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if end do @@ -470,8 +470,8 @@ contains if (dict_has_key(surface_dict, id)) then t % surface_bins(j) % scalar = dict_get_key(surface_dict, id) else - message = "Could not find surface " // trim(int_to_str(id)) // & - " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find surface " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if end do @@ -486,8 +486,8 @@ contains if (dict_has_key(universe_dict, id)) then t % universe_bins(j) % scalar = dict_get_key(universe_dict, id) else - message = "Could not find universe " // trim(int_to_str(id)) // & - " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find universe " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if end do @@ -502,8 +502,8 @@ contains if (dict_has_key(material_dict, id)) then t % material_bins(j) % scalar = dict_get_key(material_dict, id) else - message = "Could not find material " // trim(int_to_str(id)) // & - " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find material " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if end do @@ -518,8 +518,8 @@ contains if (dict_has_key(cell_dict, id)) then t % cellborn_bins(j) % scalar = dict_get_key(cell_dict, id) else - message = "Could not find material " // trim(int_to_str(id)) // & - " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find material " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if end do @@ -533,8 +533,8 @@ contains if (dict_has_key(mesh_dict, id)) then t % mesh = dict_get_key(mesh_dict, id) else - message = "Could not find mesh " // trim(int_to_str(id)) // & - " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find mesh " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if end if @@ -626,7 +626,7 @@ contains if (.not. (all(mat%atom_percent > ZERO) .or. & all(mat%atom_percent < ZERO))) then message = "Cannot mix atom and weight percents in material " // & - int_to_str(mat % id) + to_str(mat % id) call fatal_error() end if diff --git a/src/input_xml.F90 b/src/input_xml.F90 index f96682724a..e6c46cba60 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -8,7 +8,7 @@ module input_xml use global use mesh_header, only: StructuredMesh use output, only: write_message - use string, only: lower_case, int_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 @@ -224,7 +224,7 @@ contains ! Check to make sure that either material or fill was specified if (c % material == 0 .and. c % fill == 0) then message = "Neither material nor fill was specified for cell " // & - trim(int_to_str(c % id)) + trim(to_str(c % id)) call fatal_error() end if @@ -238,7 +238,7 @@ contains ! Check to make sure that surfaces were specified if (.not. associated(cell_(i) % surfaces)) then message = "No surfaces specified for cell " // & - trim(int_to_str(c % id)) + trim(to_str(c % id)) call fatal_error() end if @@ -334,11 +334,11 @@ contains n = size(surface_(i) % coeffs) if (n < coeffs_reqd) then message = "Not enough coefficients specified for surface: " // & - trim(int_to_str(s % id)) + trim(to_str(s % id)) call fatal_error() elseif (n > coeffs_reqd) then message = "Too many coefficients specified for surface: " // & - trim(int_to_str(s % id)) + trim(to_str(s % id)) call fatal_error() else allocate(s % coeffs(n)) @@ -359,7 +359,7 @@ contains s % bc = BC_PERIODIC case default message = "Unknown boundary condition '" // trim(word) // & - "' specified on surface " // trim(int_to_str(s % id)) + "' specified on surface " // trim(to_str(s % id)) call fatal_error() end select @@ -507,14 +507,14 @@ contains m % density = 1.0e-24 * val case default message = "Unkwown units '" // trim(material_(i) % density % units) & - // "' specified on material " // trim(int_to_str(m % id)) + // "' specified on material " // trim(to_str(m % id)) call fatal_error() end select ! Check to ensure material has at least one nuclide if (.not. associated(material_(i) % nuclides)) then message = "No nuclides specified on material " // & - trim(int_to_str(m % id)) + trim(to_str(m % id)) call fatal_error() end if @@ -534,7 +534,7 @@ contains ! Check for empty name on nuclide if (len_trim(nuc % name) == 0) then message = "No name specified on nuclide in material " // & - trim(int_to_str(m % id)) + trim(to_str(m % id)) call fatal_error() end if @@ -542,7 +542,7 @@ contains if (len_trim(nuc % xs) == 0) then if (default_xs == '') then message = "No cross section specified for nuclide in material " & - // trim(int_to_str(m % id)) + // trim(to_str(m % id)) call fatal_error() else nuc % xs = default_xs @@ -733,7 +733,7 @@ contains if (len_trim(tally_(i) % filters % cell) > 0 .and. & len_trim(tally_(i) % filters % surface) > 0) then message = "Cannot specify both cell and surface filters for tally " & - // trim(int_to_str(t % id)) + // trim(to_str(t % id)) call fatal_error() end if @@ -788,8 +788,8 @@ contains i_mesh = dict_get_key(mesh_dict, id) m => meshes(i_mesh) else - message = "Could not find mesh " // trim(int_to_str(id)) // & - " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find mesh " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if diff --git a/src/interpolation.F90 b/src/interpolation.F90 index 8e3f083ff5..549b720e0b 100644 --- a/src/interpolation.F90 +++ b/src/interpolation.F90 @@ -5,7 +5,7 @@ module interpolation use error, only: fatal_error use global, only: message use search, only: binary_search - use string, only: int_to_str + use string, only: to_str implicit none @@ -118,7 +118,7 @@ contains r = (log(x) - log(x0))/(log(x1) - log(x0)) y = exp((1-r)*log(y0) + r*log(y1)) case default - message = "Unsupported interpolation scheme: " // int_to_str(interp) + message = "Unsupported interpolation scheme: " // to_str(interp) call fatal_error() end select @@ -204,7 +204,7 @@ contains r = (log(x) - log(x0))/(log(x1) - log(x0)) y = exp((1-r)*log(y0) + r*log(y1)) case default - message = "Unsupported interpolation scheme: " // int_to_str(interp) + message = "Unsupported interpolation scheme: " // to_str(interp) call fatal_error() end select diff --git a/src/main.F90 b/src/main.F90 index 9586f8bc40..199da60dc0 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -10,7 +10,7 @@ program main use physics, only: transport use random_lcg, only: set_particle_seed use source, only: get_source_particle - use string, only: int_to_str + use string, only: to_str use tally, only: synchronize_tallies, write_tallies, & tally_statistics use timing, only: timer_start, timer_stop @@ -88,7 +88,7 @@ contains ! Start timer for computation call timer_start(time_compute) - message = "Simulating cycle " // trim(int_to_str(i_cycle)) // "..." + message = "Simulating cycle " // trim(to_str(i_cycle)) // "..." call write_message(8) ! Set all tallies to zero diff --git a/src/output.F90 b/src/output.F90 index b09b97fb52..7b2898528c 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -10,7 +10,7 @@ module output use global use mesh_header, only: StructuredMesh use particle_header, only: Particle, LocalCoord - use string, only: upper_case, int_to_str, real_to_str + use string, only: upper_case, to_str use tally_header, only: TallyObject implicit none @@ -67,9 +67,9 @@ contains ! Write information on number of processors #ifdef MPI write(UNIT=OUTPUT_UNIT, FMT='(1X,A)') ' MPI Processes: ' // & - trim(int_to_str(n_procs)) + trim(to_str(n_procs)) write(UNIT=UNIT_SUMMARY, FMT='(1X,"MPI Processes:",1X,A)') & - trim(int_to_str(n_procs)) + trim(to_str(n_procs)) #endif end subroutine title @@ -213,13 +213,13 @@ contains ! display type of particle select case (p % type) case (NEUTRON) - write(ou,*) 'Neutron ' // int_to_str(p % id) + write(ou,*) 'Neutron ' // to_str(p % id) case (PHOTON) - write(ou,*) 'Photon ' // int_to_str(p % id) + write(ou,*) 'Photon ' // to_str(p % id) case (ELECTRON) - write(ou,*) 'Electron ' // int_to_str(p % id) + write(ou,*) 'Electron ' // to_str(p % id) case default - write(ou,*) 'Unknown Particle ' // int_to_str(p % id) + write(ou,*) 'Unknown Particle ' // to_str(p % id) end select ! loop through each level of universes @@ -227,26 +227,26 @@ contains i = 0 do while(associated(coord)) ! Print level - write(ou,*) ' Level ' // trim(int_to_str(i)) + write(ou,*) ' Level ' // trim(to_str(i)) ! Print cell for this level if (coord % cell /= NONE) then c => cells(coord % cell) - write(ou,*) ' Cell = ' // trim(int_to_str(c % id)) + write(ou,*) ' Cell = ' // trim(to_str(c % id)) end if ! Print universe for this level if (coord % universe /= NONE) then u => universes(coord % universe) - write(ou,*) ' Universe = ' // trim(int_to_str(u % id)) + write(ou,*) ' Universe = ' // trim(to_str(u % id)) end if ! Print information on lattice if (coord % lattice /= NONE) then l => lattices(coord % lattice) - write(ou,*) ' Lattice = ' // trim(int_to_str(l % id)) - write(ou,*) ' Lattice position = (' // trim(int_to_str(& - p % coord % lattice_x)) // ',' // trim(int_to_str(& + write(ou,*) ' Lattice = ' // trim(to_str(l % id)) + write(ou,*) ' Lattice position = (' // trim(to_str(& + p % coord % lattice_x)) // ',' // trim(to_str(& p % coord % lattice_y)) // ')' end if @@ -261,13 +261,13 @@ contains ! Print surface if (p % surface /= NONE) then s => surfaces(p % surface) - write(ou,*) ' Surface = ' // int_to_str(s % id) + write(ou,*) ' Surface = ' // to_str(s % id) end if - write(ou,*) ' Weight = ' // real_to_str(p % wgt) - write(ou,*) ' Energy = ' // real_to_str(p % E) - write(ou,*) ' IE = ' // int_to_str(p % IE) - write(ou,*) ' Interpolation factor = ' // real_to_str(p % interp) + write(ou,*) ' Weight = ' // to_str(p % wgt) + write(ou,*) ' Energy = ' // to_str(p % E) + write(ou,*) ' IE = ' // to_str(p % IE) + write(ou,*) ' Interpolation factor = ' // to_str(p % interp) write(ou,*) end subroutine print_particle @@ -281,12 +281,12 @@ contains type(Reaction), pointer :: rxn write(ou,*) 'Reaction ' // reaction_name(rxn % MT) - write(ou,*) ' MT = ' // int_to_str(rxn % MT) - write(ou,*) ' Q-value = ' // real_to_str(rxn % Q_value) - write(ou,*) ' TY = ' // int_to_str(rxn % TY) - write(ou,*) ' Starting index = ' // int_to_str(rxn % IE) + write(ou,*) ' MT = ' // to_str(rxn % MT) + write(ou,*) ' Q-value = ' // to_str(rxn % Q_value) + write(ou,*) ' TY = ' // to_str(rxn % TY) + write(ou,*) ' Starting index = ' // to_str(rxn % IE) if (rxn % has_energy_dist) then - write(ou,*) ' Energy: Law ' // int_to_str(rxn % edist % law) + write(ou,*) ' Energy: Law ' // to_str(rxn % edist % law) end if write(ou,*) @@ -315,28 +315,28 @@ contains unit_ = OUTPUT_UNIT end if - write(unit_,*) 'Cell ' // int_to_str(c % id) + write(unit_,*) 'Cell ' // to_str(c % id) temp = dict_get_key(cell_dict, c % id) - write(unit_,*) ' Array Index = ' // int_to_str(temp) + write(unit_,*) ' Array Index = ' // to_str(temp) u => universes(c % universe) - write(unit_,*) ' Universe = ' // int_to_str(u % id) + write(unit_,*) ' Universe = ' // to_str(u % id) select case (c % type) case (CELL_NORMAL) write(unit_,*) ' Fill = NONE' case (CELL_FILL) u => universes(c % fill) - write(unit_,*) ' Fill = Universe ' // int_to_str(u % id) + write(unit_,*) ' Fill = Universe ' // to_str(u % id) case (CELL_LATTICE) l => lattices(c % fill) - write(unit_,*) ' Fill = Lattice ' // int_to_str(l % id) + write(unit_,*) ' Fill = Lattice ' // to_str(l % id) end select if (c % material == 0) then write(unit_,*) ' Material = NONE' else m => materials(c % material) - write(unit_,*) ' Material = ' // int_to_str(m % id) + write(unit_,*) ' Material = ' // to_str(m % id) end if - write(unit_,*) ' Parent Cell = ' // int_to_str(c % parent) + write(unit_,*) ' Parent Cell = ' // to_str(c % parent) string = "" do i = 1, c % n_surfaces select case (c % surfaces(i)) @@ -349,7 +349,7 @@ contains case (OP_DIFFERENCE) string = trim(string) // ' !' case default - string = trim(string) // ' ' // int_to_str(c % surfaces(i)) + string = trim(string) // ' ' // to_str(c % surfaces(i)) end select end do write(unit_,*) ' Surface Specification:' // trim(string) @@ -377,12 +377,12 @@ contains unit_ = OUTPUT_UNIT end if - write(unit_,*) 'Universe ' // int_to_str(univ % id) - write(unit_,*) ' Level = ' // int_to_str(univ % level) + write(unit_,*) 'Universe ' // to_str(univ % id) + write(unit_,*) ' Level = ' // to_str(univ % level) string = "" do i = 1, univ % n_cells c => cells(univ % cells(i)) - string = trim(string) // ' ' // int_to_str(c % id) + string = trim(string) // ' ' // to_str(c % id) end do write(unit_,*) ' Cells =' // trim(string) write(unit_,*) @@ -406,13 +406,13 @@ contains unit_ = OUTPUT_UNIT end if - write(unit_,*) 'Lattice ' // int_to_str(lat % id) - write(unit_,*) ' n_x = ' // int_to_str(lat % n_x) - write(unit_,*) ' n_y = ' // int_to_str(lat % n_y) - write(unit_,*) ' x0 = ' // real_to_str(lat % x0) - write(unit_,*) ' y0 = ' // real_to_str(lat % y0) - write(unit_,*) ' width_x = ' // real_to_str(lat % width_x) - write(unit_,*) ' width_y = ' // real_to_str(lat % width_y) + write(unit_,*) 'Lattice ' // to_str(lat % id) + write(unit_,*) ' n_x = ' // to_str(lat % n_x) + write(unit_,*) ' n_y = ' // to_str(lat % n_y) + write(unit_,*) ' x0 = ' // to_str(lat % x0) + write(unit_,*) ' y0 = ' // to_str(lat % y0) + write(unit_,*) ' width_x = ' // to_str(lat % width_x) + write(unit_,*) ' width_y = ' // to_str(lat % width_y) write(unit_,*) end subroutine print_lattice @@ -436,7 +436,7 @@ contains unit_ = OUTPUT_UNIT end if - write(unit_,*) 'Surface ' // int_to_str(surf % id) + write(unit_,*) 'Surface ' // to_str(surf % id) select case (surf % type) case (SURF_PX) string = "X Plane" @@ -465,14 +465,14 @@ contains string = "" do i = 1, size(surf % coeffs) - string = trim(string) // ' ' // real_to_str(surf % coeffs(i), 4) + string = trim(string) // ' ' // to_str(surf % coeffs(i), 4) end do write(unit_,*) ' Coefficients = ' // trim(string) string = "" if (allocated(surf % neighbor_pos)) then do i = 1, size(surf % neighbor_pos) - string = trim(string) // ' ' // int_to_str(surf % neighbor_pos(i)) + string = trim(string) // ' ' // to_str(surf % neighbor_pos(i)) end do end if write(unit_,*) ' Positive Neighbors = ' // trim(string) @@ -480,7 +480,7 @@ contains string = "" if (allocated(surf % neighbor_neg)) then do i = 1, size(surf % neighbor_neg) - string = trim(string) // ' ' // int_to_str(surf % neighbor_neg(i)) + string = trim(string) // ' ' // to_str(surf % neighbor_neg(i)) end do end if write(unit_,*) ' Negative Neighbors =' // trim(string) @@ -520,10 +520,10 @@ contains end if ! Write identifier for material - write(unit_,*) 'Material ' // int_to_str(mat % id) + write(unit_,*) 'Material ' // to_str(mat % id) ! Write total atom density in atom/b-cm - write(unit_,*) ' Atom Density = ' // trim(real_to_str(mat % density)) & + write(unit_,*) ' Atom Density = ' // trim(to_str(mat % density)) & // ' atom/b-cm' ! Write atom density for each nuclide in material @@ -532,7 +532,7 @@ contains nuc => nuclides(mat % nuclide(i)) density = mat % atom_density(i) string = ' ' // trim(nuc % name) // ' = ' // & - trim(real_to_str(density)) // ' atom/b-cm' + trim(to_str(density)) // ' atom/b-cm' write(unit_,*) trim(string) end do @@ -569,14 +569,14 @@ contains unit_ = OUTPUT_UNIT end if - write(unit_,*) 'Tally ' // int_to_str(t % id) + write(unit_,*) 'Tally ' // to_str(t % id) if (t % n_bins(T_CELL) > 0) then string = "" do i = 1, t % n_bins(T_CELL) id = t % cell_bins(i) % scalar c => cells(id) - string = trim(string) // ' ' // trim(int_to_str(c % id)) + string = trim(string) // ' ' // trim(to_str(c % id)) end do write(unit_, *) ' Cell Bins:' // trim(string) end if @@ -586,7 +586,7 @@ contains do i = 1, t % n_bins(T_SURFACE) id = t % surface_bins(i) % scalar s => surfaces(id) - string = trim(string) // ' ' // trim(int_to_str(s % id)) + string = trim(string) // ' ' // trim(to_str(s % id)) end do write(unit_, *) ' Surface Bins:' // trim(string) end if @@ -596,7 +596,7 @@ contains do i = 1, t % n_bins(T_UNIVERSE) id = t % universe_bins(i) % scalar u => universes(id) - string = trim(string) // ' ' // trim(int_to_str(u % id)) + string = trim(string) // ' ' // trim(to_str(u % id)) end do write(unit_, *) ' Material Bins:' // trim(string) end if @@ -606,7 +606,7 @@ contains do i = 1, t % n_bins(T_MATERIAL) id = t % material_bins(i) % scalar m => materials(id) - string = trim(string) // ' ' // trim(int_to_str(m % id)) + string = trim(string) // ' ' // trim(to_str(m % id)) end do write(unit_, *) ' Material Bins:' // trim(string) end if @@ -615,9 +615,9 @@ contains string = "" id = t % mesh sm => meshes(id) - string = trim(string) // ' ' // trim(int_to_str(sm % dimension(1))) + string = trim(string) // ' ' // trim(to_str(sm % dimension(1))) do i = 2, sm % n_dimension - string = trim(string) // ' x ' // trim(int_to_str(sm % dimension(i))) + string = trim(string) // ' x ' // trim(to_str(sm % dimension(i))) end do write(unit_, *) ' Mesh Bins:' // trim(string) end if @@ -627,7 +627,7 @@ contains do i = 1, t % n_bins(T_CELLBORN) id = t % cellborn_bins(i) % scalar c => cells(id) - string = trim(string) // ' ' // trim(int_to_str(c % id)) + string = trim(string) // ' ' // trim(to_str(c % id)) end do write(unit_, *) ' Birth Region Bins:' // trim(string) end if @@ -635,7 +635,7 @@ contains if (t % n_bins(T_ENERGYIN) > 0) then string = "" do i = 1, t % n_bins(T_ENERGYIN) + 1 - string = trim(string) // ' ' // trim(real_to_str(& + string = trim(string) // ' ' // trim(to_str(& t % energy_in(i))) end do write(unit_,*) ' Incoming Energy Bins:' // trim(string) @@ -644,7 +644,7 @@ contains if (t % n_bins(T_ENERGYOUT) > 0) then string = "" do i = 1, t % n_bins(T_ENERGYOUT) + 1 - string = trim(string) // ' ' // trim(real_to_str(& + string = trim(string) // ' ' // trim(to_str(& t % energy_out(i))) end do write(unit_,*) ' Outgoing Energy Bins:' // trim(string) @@ -751,14 +751,14 @@ contains ! Basic nuclide information write(unit_,*) 'Nuclide ' // trim(nuc % name) - write(unit_,*) ' zaid = ' // trim(int_to_str(nuc % zaid)) - write(unit_,*) ' awr = ' // trim(real_to_str(nuc % awr)) - write(unit_,*) ' kT = ' // trim(real_to_str(nuc % kT)) - write(unit_,*) ' # of grid points = ' // trim(int_to_str(nuc % n_grid)) + write(unit_,*) ' zaid = ' // trim(to_str(nuc % zaid)) + write(unit_,*) ' awr = ' // trim(to_str(nuc % awr)) + write(unit_,*) ' kT = ' // trim(to_str(nuc % kT)) + write(unit_,*) ' # of grid points = ' // trim(to_str(nuc % n_grid)) write(unit_,*) ' Fissionable = ', nuc % fissionable - write(unit_,*) ' # of fission reactions = ' // trim(int_to_str(nuc % n_fission)) - write(unit_,*) ' # of reactions = ' // trim(int_to_str(nuc % n_reaction)) - write(unit_,*) ' Size of cross sections = ' // trim(int_to_str(& + write(unit_,*) ' # of fission reactions = ' // trim(to_str(nuc % n_fission)) + write(unit_,*) ' # of reactions = ' // trim(to_str(nuc % n_reaction)) + write(unit_,*) ' Size of cross sections = ' // trim(to_str(& size_xs)) // ' bytes' write(unit_,*) ' Reaction Q-value Mult IE size(angle) size(energy)' @@ -792,18 +792,18 @@ contains if (nuc % urr_present) then urr => nuc % urr_data write(unit_,*) ' Unresolved resonance probability table:' - write(unit_,*) ' # of energies = ' // trim(int_to_str(urr % n_energy)) - write(unit_,*) ' # of probabilities = ' // trim(int_to_str(urr % n_prob)) - write(unit_,*) ' Interpolation = ' // trim(int_to_str(urr % interp)) - write(unit_,*) ' Inelastic flag = ' // trim(int_to_str(urr % inelastic_flag)) - write(unit_,*) ' Absorption flag = ' // trim(int_to_str(urr % absorption_flag)) + write(unit_,*) ' # of energies = ' // trim(to_str(urr % n_energy)) + write(unit_,*) ' # of probabilities = ' // trim(to_str(urr % n_prob)) + write(unit_,*) ' Interpolation = ' // trim(to_str(urr % interp)) + write(unit_,*) ' Inelastic flag = ' // trim(to_str(urr % inelastic_flag)) + write(unit_,*) ' Absorption flag = ' // trim(to_str(urr % absorption_flag)) write(unit_,*) ' Multiply by smooth? ', urr % multiply_smooth - write(unit_,*) ' Min energy = ', trim(real_to_str(urr % energy(1))) - write(unit_,*) ' Max energy = ', trim(real_to_str(urr % energy(urr % n_energy))) + write(unit_,*) ' Min energy = ', trim(to_str(urr % energy(1))) + write(unit_,*) ' Max energy = ', trim(to_str(urr % energy(urr % n_energy))) end if ! Write total memory used - write(unit_,*) ' Total memory used = ' // trim(int_to_str(size_total)) & + write(unit_,*) ' Total memory used = ' // trim(to_str(size_total)) & // ' bytes' ! Blank line at end of nuclide @@ -861,8 +861,8 @@ contains ! print summary of unionized energy grid call header("UNIONIZED ENERGY GRID", unit=UNIT_SUMMARY) - write(UNIT_SUMMARY,*) "Points on energy grid: " // trim(int_to_str(n_grid)) - write(UNIT_SUMMARY,*) "Extra storage required: " // trim(int_to_str(& + write(UNIT_SUMMARY,*) "Points on energy grid: " // trim(to_str(n_grid)) + write(UNIT_SUMMARY,*) "Extra storage required: " // trim(to_str(& n_grid*n_nuclides_total*4)) // " bytes" ! print summary of variance reduction @@ -872,9 +872,9 @@ contains else write(UNIT_SUMMARY,100) "Survival Biasing:", "off" end if - string = real_to_str(weight_cutoff) + string = to_str(weight_cutoff) write(UNIT_SUMMARY,100) "Weight Cutoff:", trim(string) - string = real_to_str(weight_survive) + string = to_str(weight_survive) write(UNIT_SUMMARY,100) "Survival weight:", trim(string) ! Format descriptor for columns @@ -893,16 +893,16 @@ contains call header("PLOTTING SUMMARY") ! Print plotting origin - write(ou,100) "Plotting Origin:", trim(real_to_str(plot_origin(1))) // & - " " // trim(real_to_str(plot_origin(2))) // " " // & - trim(real_to_str(plot_origin(3))) + write(ou,100) "Plotting Origin:", trim(to_str(plot_origin(1))) // & + " " // trim(to_str(plot_origin(2))) // " " // & + trim(to_str(plot_origin(3))) ! Print plotting width - write(ou,100) "Plotting Width:", trim(real_to_str(plot_width(1))) // & - " " // trim(real_to_str(plot_width(2))) + write(ou,100) "Plotting Width:", trim(to_str(plot_width(1))) // & + " " // trim(to_str(plot_width(2))) ! Print pixel width - write(ou,100) "Pixel Width:", trim(real_to_str(pixel)) + write(ou,100) "Pixel Width:", trim(to_str(pixel)) write(ou,*) ! Format descriptor for columns @@ -944,7 +944,7 @@ contains ! display calculate rate and final keff total_particles = n_particles * n_cycles speed = real(total_particles) / time_compute % elapsed - string = real_to_str(speed) + string = to_str(speed) write(ou,101) "Calculation Rate", trim(string) write(ou,102) "Final Keff", keff, keff_std write(ou,*) diff --git a/src/physics.F90 b/src/physics.F90 index 449ee990b5..685e379bd4 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -15,7 +15,7 @@ module physics use particle_header, only: Particle, LocalCoord use random_lcg, only: prn use search, only: binary_search - use string, only: int_to_str, real_to_str + use string, only: to_str use tally, only: score_tally, score_surface_current implicit none @@ -54,12 +54,12 @@ contains end if if (verbosity >= 9 .or. trace) then - message = "Simulating Particle " // trim(int_to_str(p % id)) + message = "Simulating Particle " // trim(to_str(p % id)) call write_message() end if if (verbosity >= 10 .or. trace) then - message = " Born in cell " // trim(int_to_str(& + message = " Born in cell " // trim(to_str(& cells(p % coord % cell) % id)) call write_message() end if @@ -126,7 +126,7 @@ contains ! If particle has too many events, display warning and kill it n_event = n_event + 1 if (n_event == MAX_EVENTS) then - message = "Particle " // trim(int_to_str(p%id)) // " underwent " & + message = "Particle " // trim(to_str(p%id)) // " underwent " & // "maximum number of events." call warning() p % alive = .false. @@ -168,7 +168,7 @@ contains ! Display information about collision if (verbosity >= 10 .or. trace) then message = " " // trim(reaction_name(MT)) // ". Energy = " // & - trim(real_to_str(p % E * 1e6_8)) // " eV." + trim(to_str(p % E * 1e6_8)) // " eV." call write_message() end if @@ -417,7 +417,7 @@ contains if (i > nuc % n_reaction) then message = "Did not sample any reaction for nuclide " // & trim(nuc % name) // " on material " // & - trim(int_to_str(mat % id)) + trim(to_str(mat % id)) call fatal_error() end if @@ -1136,7 +1136,7 @@ contains mu = mu0 + (sqrt(p0*p0 + 2*frac*(xi - c_k))-p0)/frac end if else - message = "Unknown interpolation type: " // trim(int_to_str(interp)) + message = "Unknown interpolation type: " // trim(to_str(interp)) call fatal_error() end if @@ -1148,7 +1148,7 @@ contains end if else - message = "Unknown angular distribution type: " // trim(int_to_str(type)) + message = "Unknown angular distribution type: " // trim(to_str(type)) call fatal_error() end if @@ -1455,7 +1455,7 @@ contains p_l_k)/frac end if else - message = "Unknown interpolation type: " // trim(int_to_str(INTT)) + message = "Unknown interpolation type: " // trim(to_str(INTT)) call fatal_error() end if @@ -1696,7 +1696,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 - message = "Unknown interpolation type: " // trim(int_to_str(INTT)) + message = "Unknown interpolation type: " // trim(to_str(INTT)) call fatal_error() end if @@ -1826,7 +1826,7 @@ contains p_l_k)/frac end if else - message = "Unknown interpolation type: " // trim(int_to_str(INTT)) + message = "Unknown interpolation type: " // trim(to_str(INTT)) call fatal_error() end if @@ -1879,7 +1879,7 @@ contains mu_out = mu_k + (sqrt(p_k*p_k + 2*frac*(r3 - c_k))-p_k)/frac end if else - message = "Unknown interpolation type: " // trim(int_to_str(JJ)) + message = "Unknown interpolation type: " // trim(to_str(JJ)) call fatal_error() end if diff --git a/src/source.F90 b/src/source.F90 index f296c2c5b6..1ecea63eda 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -9,7 +9,7 @@ module source use particle_header, only: Particle, initialize_particle use physics, only: watt_spectrum use random_lcg, only: prn, set_particle_seed - use string, only: int_to_str + use string, only: to_str implicit none @@ -51,7 +51,7 @@ contains bytes = maxwork * 64 / 8 #endif message = "Could not allocate source bank. Attempted to allocate " & - // trim(int_to_str(bytes)) // " bytes." + // trim(to_str(bytes)) // " bytes." call fatal_error() end if @@ -64,7 +64,7 @@ contains bytes = 3 * maxwork * 64 / 8 #endif message = "Could not allocate fission bank. Attempted to allocate " & - // trim(int_to_str(bytes)) // " bytes." + // trim(to_str(bytes)) // " bytes." call fatal_error() end if @@ -74,7 +74,7 @@ contains p_max = external_source%values(4:6) else message = "Unsupported external source type: " // & - int_to_str(external_source%type) + to_str(external_source%type) call fatal_error() end if diff --git a/src/string.F90 b/src/string.F90 index ae8991b849..2800d84c12 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -6,8 +6,8 @@ module string implicit none - interface int_to_str - module procedure int4_to_str, int8_to_str + interface to_str + module procedure int4_to_str, int8_to_str, real_to_str end interface contains diff --git a/src/tally.F90 b/src/tally.F90 index 23d0e9bc9e..a1d66674ec 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -7,7 +7,7 @@ module tally use mesh_header, only: StructuredMesh use output, only: header use search, only: binary_search - use string, only: int_to_str, real_to_str + use string, only: to_str use tally_header, only: TallyScore, TallyMapItem, TallyMapElement #ifdef MPI @@ -966,7 +966,7 @@ contains t => tallies(i) ! Write header block - call header("TALLY " // trim(int_to_str(t % id)), unit=UNIT_TALLY, level=3) + call header("TALLY " // trim(to_str(t % id)), unit=UNIT_TALLY, level=3) ! Handle surface current tallies separately if (t % surface_current) then @@ -1048,8 +1048,8 @@ contains do k = 1, t % n_macro_bins write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), macro_name(abs(t % macro_bins(k) % scalar)), & - real_to_str(t % scores(score_index,k) % val), & - trim(real_to_str(t % scores(score_index,k) % val_sq)) + to_str(t % scores(score_index,k) % val), & + trim(to_str(t % scores(score_index,k) % val_sq)) end do indent = indent - 2 @@ -1099,14 +1099,14 @@ contains end if do i = 1, m % dimension(1) - string = "Mesh Index (" // trim(int_to_str(i)) // ", " + string = "Mesh Index (" // trim(to_str(i)) // ", " len1 = len_trim(string) do j = 1, m % dimension(2) - string = string(1:len1+1) // trim(int_to_str(j)) // ", " + string = string(1:len1+1) // trim(to_str(j)) // ", " len2 = len_trim(string) do k = 1, m % dimension(3) ! Write mesh cell index - string = string(1:len2+1) // trim(int_to_str(k)) // ")" + string = string(1:len2+1) // trim(to_str(k)) // ")" write(UNIT=UNIT_TALLY, FMT='(1X,A)') trim(string) do l = 1, n @@ -1125,15 +1125,15 @@ contains score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Left", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) bins(TS_SURFACE) = OUT_RIGHT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Left", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) ! Right Surface bins(1:3) = (/ i, j, k /) + 1 @@ -1141,15 +1141,15 @@ contains score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Right", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) bins(TS_SURFACE) = OUT_RIGHT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Right", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) ! Back Surface bins(1:3) = (/ i, j-1, k /) + 1 @@ -1157,15 +1157,15 @@ contains score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Back", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) bins(TS_SURFACE) = OUT_FRONT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Back", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) ! Front Surface bins(1:3) = (/ i, j, k /) + 1 @@ -1173,15 +1173,15 @@ contains score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Front", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) bins(TS_SURFACE) = OUT_FRONT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Front", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) ! Bottom Surface bins(1:3) = (/ i, j, k-1 /) + 1 @@ -1189,15 +1189,15 @@ contains score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Bottom", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) bins(TS_SURFACE) = OUT_TOP score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Bottom", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) ! Top Surface bins(1:3) = (/ i, j, k /) + 1 @@ -1205,15 +1205,15 @@ contains score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Top", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) bins(TS_SURFACE) = OUT_TOP score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Top", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) end do end do @@ -1243,38 +1243,38 @@ contains select case(filter_type) case (T_UNIVERSE) i = t % universe_bins(bin) % scalar - label = int_to_str(universes(i) % id) + label = to_str(universes(i) % id) case (T_MATERIAL) i = t % material_bins(bin) % scalar - label = int_to_str(materials(i) % id) + label = to_str(materials(i) % id) case (T_CELL) i = t % cell_bins(bin) % scalar - label = int_to_str(cells(i) % id) + label = to_str(cells(i) % id) case (T_CELLBORN) i = t % cellborn_bins(bin) % scalar - label = int_to_str(cells(i) % id) + label = to_str(cells(i) % id) case (T_SURFACE) i = t % surface_bins(bin) % scalar - label = int_to_str(surfaces(i) % id) + label = to_str(surfaces(i) % id) case (T_MESH) m => meshes(t % mesh) allocate(ijk(m % n_dimension)) call bin_to_mesh_indices(m, bin, ijk) if (m % n_dimension == 2) then - label = "Index (" // trim(int_to_str(ijk(1))) // ", " // & - trim(int_to_str(ijk(2))) // ")" + label = "Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ")" elseif (m % n_dimension == 3) then - label = "Index (" // trim(int_to_str(ijk(1))) // ", " // & - trim(int_to_str(ijk(2))) // ", " // trim(int_to_str(ijk(3))) // ")" + label = "Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" end if case (T_ENERGYIN) E0 = t % energy_in(bin) E1 = t % energy_in(bin + 1) - label = "[" // trim(real_to_str(E0)) // ", " // trim(real_to_str(E1)) // ")" + label = "[" // trim(to_str(E0)) // ", " // trim(to_str(E1)) // ")" case (T_ENERGYOUT) E0 = t % energy_out(bin) E1 = t % energy_out(bin + 1) - label = "[" // trim(real_to_str(E0)) // ", " // trim(real_to_str(E1)) // ")" + label = "[" // trim(to_str(E0)) // ", " // trim(to_str(E1)) // ")" end select end function get_label From 5964082dde208d71c1d744daf1f9dcfbf50191a1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 21 Dec 2011 00:48:46 -0500 Subject: [PATCH 32/33] Added particle id to fatal error messages. --- src/geometry.F90 | 21 +++++++++++---------- src/physics.F90 | 3 +-- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 983cf5d3c1..35612f23cc 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -285,7 +285,8 @@ contains ! Do not handle reflective boundary conditions on lower universes if (p % in_lower_universe) then - message = "Cannot reflect particle off surface in a lower universe." + message = "Cannot reflect particle " // trim(to_str(p % id)) // & + " off surface in a lower universe." call fatal_error() end if @@ -428,9 +429,9 @@ contains ! Couldn't find next cell anywhere! if ((.not. found) .and. (.not. plotting)) then - message = "After particle crossed surface " // trim(to_str( & - surfaces(abs(p%surface)) % id)) // " it could not be located in & - &any cell and it did not leak." + message = "After particle " // trim(to_str(p % id)) // " crossed surface " & + // trim(to_str(surfaces(abs(p%surface)) % id)) // " it could not be & + &located in any cell and it did not leak." call fatal_error() end if @@ -495,14 +496,14 @@ contains i_x = p % coord % lattice_x i_y = p % coord % lattice_y if (i_x < 1 .or. i_x > lat % n_x) then - message = "Reached edge of lattice " // trim(to_str(lat % id)) // & - " at position (" // trim(to_str(i_x)) // "," // & - trim(to_str(i_y)) // ")." + 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 = "Reached edge of lattice " // trim(to_str(lat % id)) // & - " at position (" // trim(to_str(i_x)) // "," // & - trim(to_str(i_y)) // ")." + 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 diff --git a/src/physics.F90 b/src/physics.F90 index 685e379bd4..d3bf1fc894 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -44,8 +44,7 @@ contains call find_cell(p, found_cell) ! Particle couldn't be located if (.not. found_cell) then - write(message, '(A,3ES11.3)') & - "Could not locate cell for particle at: ", p % coord0 % xyz + message = "Could not locate particle " // trim(to_str(p % id)) call fatal_error() end if From 44ebfe0a790e808eec76e83263ce17e81825ba4e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 21 Dec 2011 11:06:47 -0500 Subject: [PATCH 33/33] Fixed bug in selection of energy for probability tables. --- src/cross_section.F90 | 2 +- src/physics.F90 | 20 ++++++++++---------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 79ac18d00c..367d277cf5 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -296,7 +296,7 @@ contains ! determine energy table i_energy = 1 do - if (p % E > urr % energy(i_energy)) exit + if (p % E < urr % energy(i_energy + 1)) exit i_energy = i_energy + 1 end do diff --git a/src/physics.F90 b/src/physics.F90 index d3bf1fc894..869acc367a 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -30,14 +30,14 @@ contains type(Particle), pointer :: p - integer :: surface_crossed ! surface which particle is on - integer :: lattice_crossed ! lattice boundary which particle crossed - integer :: last_cell ! most recent cell particle was in - integer :: n_event ! number of collisions/crossings - real(8) :: d_boundary ! distance to nearest boundary - real(8) :: d_collision ! sampled distance to collision - real(8) :: distance ! distance particle travels - logical :: found_cell ! found cell which particle is in? + integer :: surface_crossed ! surface which particle is on + integer :: lattice_crossed ! lattice boundary which particle crossed + integer :: last_cell ! most recent cell particle was in + integer :: n_event ! number of collisions/crossings + real(8) :: d_boundary ! distance to nearest boundary + real(8) :: d_collision ! sampled distance to collision + real(8) :: distance ! distance particle travels + logical :: found_cell ! found cell which particle is in? type(LocalCoord), pointer :: coord => null() if (p % coord % cell == NONE) then @@ -125,8 +125,8 @@ contains ! If particle has too many events, display warning and kill it n_event = n_event + 1 if (n_event == MAX_EVENTS) then - message = "Particle " // trim(to_str(p%id)) // " underwent " & - // "maximum number of events." + message = "Particle " // trim(to_str(p%id)) // " underwent maximum & + &number of events." call warning() p % alive = .false. end if