From 9a60b43946d1bfb42ad5b8086a2958dda1444023 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 21 Apr 2011 20:24:08 +0000 Subject: [PATCH] Fixed bug in sampling of angle. Changed Bank type to have uid as integer(8) to allow for more than 2^32 - 1 particles to be simulated. --- ChangeLog | 18 +++++++++++++ src/geometry.f90 | 6 ++--- src/global.f90 | 45 ++++++++++++++++++++++---------- src/lattice_sample | 42 ++++++++++++++++++------------ src/main.f90 | 12 ++++++--- src/mpi_routines.f90 | 62 ++++++++++++++++++++++++++++---------------- src/output.f90 | 6 ++--- src/physics.f90 | 19 +++++++++----- src/source.f90 | 20 +++++++------- src/tallies.f90 | 4 +-- src/types.f90 | 42 +++++++++++++++--------------- 11 files changed, 177 insertions(+), 99 deletions(-) diff --git a/ChangeLog b/ChangeLog index 11501105a1..8874be2888 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,21 @@ +2011-04-21 Paul Romano + + * geometry.f90: Removed level on calls to message(). + * global.f90: Adjustments to make necessary integers into 8-byte + integers. Changed PI to non-compiler dependent constant. + * lattice_sample: Changed lattice to 20 x 20. + * main.f90: New timing data from fission bank synchronization + * mpi_routines.f90: Switching to 8-byte integers, more comments, + changed timing variables. + * output.f90: Made the level argument in message subroutine + optional. + * physics.f90: Fixed bug in sample_angle whereby mu could possibly + be greater than one (very rarely). This caused rotate_angle to + fail when given mu>1. Removed level from message calls. + * source.f90: Change necessary variables to integer(8) + * tallies.f90: Change MPI_REDUCE to use MPI_INTEGER8 + * types.f90: Change Bank type to use integer(8) for uid. + 2011-04-13 Paul Romano * main.f90: Simple timing for total time elapsed and bank diff --git a/src/geometry.f90 b/src/geometry.f90 index b7140a8200..a1f9e66831 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -185,7 +185,7 @@ contains surf => surfaces(abs(p%surface)) if (verbosity >= 10) then msg = " Crossing surface " // trim(int_to_str(surf%uid)) - call message(msg, 10) + call message(msg) end if ! check for leakage @@ -193,7 +193,7 @@ contains p%alive = .false. if (verbosity >= 10) then msg = " Leaked out of surface " // trim(int_to_str(surf%uid)) - call message(msg, 10) + call message(msg) end if return end if @@ -294,7 +294,7 @@ contains if (verbosity >= 10) then msg = " Crossing lattice" - call message(msg, 10) + call message(msg) end if lat => lattices(p % lattice) diff --git a/src/global.f90 b/src/global.f90 index b21745fd8c..7301d798af 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -50,9 +50,9 @@ module global real(8), allocatable :: e_grid(:) ! energies on unionized grid ! Histories/cycles/etc for both external source and criticality - integer :: n_particles ! # of particles (per cycle for criticality) - integer :: n_cycles ! # of cycles - integer :: n_inactive ! # of inactive cycles + integer(8) :: n_particles ! # of particles (per cycle for criticality) + integer :: n_cycles ! # of cycles + integer :: n_inactive ! # of inactive cycles ! External source type(ExtSource), target :: external_source @@ -60,10 +60,10 @@ module global ! Source and fission bank type(Particle), allocatable, target :: source_bank(:) type(Bank), allocatable, target :: fission_bank(:) - integer :: n_bank ! # of sites in fission bank - integer :: bank_first ! index of first particle in bank - integer :: bank_last ! index of last particle in bank - integer :: work ! number of particles per processor + integer(8) :: n_bank ! # of sites in fission bank + integer(8) :: bank_first ! index of first particle in bank + integer(8) :: bank_last ! index of last particle in bank + integer :: work ! number of particles per processor ! cycle keff real(8) :: keff @@ -81,7 +81,7 @@ module global ! Physical constants real(8), parameter :: & - & PI = 2.0_8*acos(0.0_8), & ! pi + & PI = 3.1415926535898_8, & ! pi & MASS_NEUTRON = 1.0086649156, & ! mass of a neutron & MASS_PROTON = 1.00727646677, & ! mass of a proton & AMU = 1.66053873e-27, & ! 1 amu in kg @@ -246,6 +246,10 @@ module global integer, parameter :: VERSION_MINOR = 2 integer, parameter :: VERSION_RELEASE = 2 + interface int_to_str + module procedure int4_to_str, int8_to_str + end interface + contains !===================================================================== @@ -300,19 +304,32 @@ contains end subroutine free_memory !===================================================================== -! INT_TO_STR converts an integer to a string. Right now, it is limited -! to integers less than 10 billion. +! INT4_TO_STR converts an integer(4) to a string. !===================================================================== - function int_to_str(num) result(str) + function int4_to_str(num) result(str) integer, intent(in) :: num - character(10) :: str + character(11) :: str - write (str, '(I10)') num + write (str, '(I11)') num str = adjustl(str) - end function int_to_str + end function int4_to_str + +!===================================================================== +! INT8_TO_STR converts an integer(8) to a string. +!===================================================================== + + function int8_to_str(num) result(str) + + integer(8), intent(in) :: num + character(21) :: str + + write (str, '(I21)') num + str = adjustl(str) + + end function int8_to_str !===================================================================== ! STR_TO_INT converts a string to an integer. diff --git a/src/lattice_sample b/src/lattice_sample index cc34d57294..8dc13c23c7 100644 --- a/src/lattice_sample +++ b/src/lattice_sample @@ -4,22 +4,32 @@ cell 1 0 fill 111 1 -2 3 -4 cell 100 3 40 -5 cell 101 3 41 5 -lattice 111 rect 10 10 -5.0 -5.0 1.0 1.0 & - 3 3 3 3 3 3 3 3 3 3 & - 3 3 3 3 3 3 3 3 3 3 & - 3 3 3 3 3 3 3 3 3 3 & - 3 3 3 3 3 3 3 3 3 3 & - 3 3 3 3 3 3 3 3 3 3 & - 3 3 3 3 3 3 3 3 3 3 & - 3 3 3 3 3 3 3 3 3 3 & - 3 3 3 3 3 3 3 3 3 3 & - 3 3 3 3 3 3 3 3 3 3 & - 3 3 3 3 3 3 3 3 3 3 +lattice 111 rect 20 20 -10.0 -10.0 1.0 1.0 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 & + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 -surface 1 px -4.9 -surface 2 px 4.9 -surface 3 py -4.9 -surface 4 py 4.9 +surface 1 px -9.9 +surface 2 px 9.9 +surface 3 py -9.9 +surface 4 py 9.9 surface 5 cylz 0 0 0.4 bc 1 vacuum @@ -37,5 +47,5 @@ material 41 -1.0 & source box -4 -4 -4 4 4 4 xs_data /opt/serpent/xsdata/endfb7 -criticality 150 25 10000 +criticality 20 10 100000 verbosity 7 diff --git a/src/main.f90 b/src/main.f90 index 025e50ee0c..5067c18c02 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -12,7 +12,7 @@ program main use cross_section, only: read_xsdata, material_total_xs use ace, only: read_xs use energy_grid, only: unionized_grid, original_indices - use mpi_routines, only: setup_mpi, synchronize_bank, t_assemble + use mpi_routines, only: setup_mpi, synchronize_bank, t_sync use tallies, only: calculate_keff #ifdef MPI @@ -116,6 +116,9 @@ contains #endif CYCLE_LOOP: do i_cycle = 1, n_cycles + + msg = "Simulating cycle " // trim(int_to_str(i_cycle)) // "..." + call message(msg, 8) ! Set all tallies to zero n_bank = 0 @@ -154,8 +157,11 @@ contains ! print run time t1 = MPI_WTIME() if (master) then - write(6,"(A,F8.2,A)") "Time elapsed = ", t1 - t0, " s" - write(6,"(A,G10.4)") "Assemble time = ", t_assemble + print *, "Time elapsed = " // real_to_str(t1 - t0) + print *, "Init time = " // real_to_str(t_sync(1)) + print *, "Sample time = " // real_to_str(t_sync(2)) + print *, "Send/recv time = " // real_to_str(t_sync(3)) + print *, "Rebuild time = " // real_to_str(t_sync(4)) end if #endif diff --git a/src/mpi_routines.f90 b/src/mpi_routines.f90 index caaa2dc67f..4c1145ded2 100644 --- a/src/mpi_routines.f90 +++ b/src/mpi_routines.f90 @@ -14,7 +14,7 @@ module mpi_routines integer :: MPI_BANK ! MPI datatype for fission bank integer(8) :: bank_index ! Fission bank site unique identifier - real(8) :: t_assemble + real(8) :: t_sync(4) contains @@ -80,11 +80,12 @@ contains ! Define MPI_BANK for fission sites bank_blocks = (/ 1, 3, 3, 1 /) - bank_types = (/ MPI_INTEGER, MPI_REAL8, MPI_REAL8, MPI_REAL8 /) + bank_types = (/ MPI_INTEGER8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /) call MPI_TYPE_CREATE_STRUCT(4, bank_blocks, bank_disp, & & bank_types, MPI_BANK, ierr) call MPI_TYPE_COMMIT(MPI_BANK, ierr) + t_sync = ZERO #else ! if no MPI, set processor to master mpi_enabled = .false. @@ -107,25 +108,24 @@ contains #ifdef MPI integer :: status(MPI_STATUS_SIZE) ! message status #endif - integer :: ierr - integer :: start ! starting index in local fission bank - integer :: finish ! ending index in local fission bank - integer :: total ! total sites in global fission bank - integer :: count ! index for source bank - integer :: index ! index for uid -- accounts for all nodes - 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 :: sites_needed ! # of sites to be sampled - integer :: sites_remaining ! # of sites left in fission bank - real(8) :: t0 - real(8) :: t1 + integer :: ierr + 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) :: count ! index for source bank + integer(8) :: index ! index for uid -- accounts for all nodes + 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 + integer(8) :: sites_remaining ! # of sites left in fission bank + real(8) :: t0, t1, t2, t3, t4 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 character(250) :: msg - msg = "Synchronizing fission bank..." + msg = "Collecting number of fission sites..." call message(msg, 8) #ifdef MPI @@ -135,11 +135,11 @@ contains ! Determine starting index for fission bank and total sites in ! fission bank start = 0 - call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER, MPI_SUM, & + call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, & & MPI_COMM_WORLD, ierr) finish = start + n_bank total = finish - call MPI_BCAST(total, 1, MPI_INTEGER, n_procs - 1, & + call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, & & MPI_COMM_WORLD, ierr) ! Check if there are no fission sites @@ -147,15 +147,18 @@ contains msg = "No fission sites banked!" call error(msg) end if + + t1 = MPI_WTIME() + t_sync(1) = t_sync(1) + (t1 - t0) #else - start = 0 + start = 0_8 finish = n_bank total = n_bank #endif allocate(temp_sites(2*work)) - count = 0 ! Index for source_bank - index = 0 ! Index for uid -- must account for all nodes + count = 0_8 ! Index for local source_bank + index = 0_8 ! Index for global source uid -- must account for all nodes if (total < n_particles) then sites_needed = mod(n_particles,total) @@ -163,6 +166,9 @@ contains sites_needed = n_particles end if + msg = "Sampling fission sites..." + call message(msg, 8) + ! ================================================================ ! SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES do i = 1, total @@ -203,6 +209,12 @@ contains send_to_right = temp_sites(1)%uid + (count-1) - bank_last #ifdef MPI + t2 = MPI_WTIME() + t_sync(2) = t_sync(2) + (t2 - t1) + + msg = "Sending fission sites..." + call message(msg, 8) + ! ================================================================ ! SEND BANK SITES TO NEIGHBORS allocate(left_bank(abs(send_to_left))) @@ -224,8 +236,14 @@ contains call MPI_SEND(temp_sites(1), send_to_left, MPI_BANK, rank-1, 1, & & MPI_COMM_WORLD, ierr) end if + + t3 = MPI_WTIME() + t_sync(3) = t_sync(3) + (t3 - t2) #endif + msg = "Constructing source bank..." + call message(msg, 8) + ! ================================================================ ! RECONSTRUCT SOURCE BANK if (send_to_left < 0 .and. send_to_right >= 0) then @@ -254,8 +272,8 @@ contains source_index = 0 #ifdef MPI - t1 = MPI_WTIME() - t_assemble = t_assemble + (t1 - t0) + t4 = MPI_WTIME() + t_sync(4) = t_sync(4) + (t4 - t3) deallocate(left_bank ) deallocate(right_bank) diff --git a/src/output.f90 b/src/output.f90 index 1ebcc786af..3fb3507979 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -98,15 +98,15 @@ contains subroutine message(msg, level) character(*), intent(in) :: msg - integer, intent(in) :: level + integer, optional :: level integer :: n_lines integer :: i ! Only allow master to print to screen - if (.not. master) return + if (.not. master .and. present(level)) return - if (level <= verbosity) then + if (.not. present(level) .or. level <= verbosity) then n_lines = (len_trim(msg)-1)/79 + 1 do i = 1, n_lines write(ou, fmt='(1X,A79)') msg(79*(i-1)+1:79*i) diff --git a/src/physics.f90 b/src/physics.f90 index 107188a35d..aa3edca1cb 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -48,12 +48,12 @@ contains if (verbosity >= 9) then msg = "Simulating Particle " // trim(int_to_str(p % uid)) - call message(msg, 9) + call message(msg) end if if (verbosity >= 10) then msg = " Born in cell " // trim(int_to_str(cells(p%cell)%uid)) - call message(msg, 10) + call message(msg) end if ! find energy index, interpolation factor @@ -224,10 +224,10 @@ contains if (r1 < prob) exit end do - if (verbosity >= 10 .or. p % uid == 4593) then + if (verbosity >= 10) then msg = " " // trim(reaction_name(rxn%MT)) // " with nuclide " // & & trim(table%name) - call message(msg, 8) + call message(msg) end if ! call appropriate subroutine @@ -783,7 +783,7 @@ contains if (verbosity >= 10) then cell_num = cells(p % cell)%uid msg = " Absorbed in cell " // trim(int_to_str(cell_num)) - call message(msg, 10) + call message(msg) end if end subroutine n_absorption @@ -897,7 +897,14 @@ contains msg = "Unknown interpolation type: " // trim(int_to_str(interp)) call error(msg) end if - + + if (abs(mu) > ONE) then + msg = "Sampled cosine of angle outside [-1, 1)." + call warning(msg) + + mu = sign(ONE,mu) + end if + else msg = "Unknown angular distribution type: " // trim(int_to_str(type)) call error(msg) diff --git a/src/source.f90 b/src/source.f90 index 6c64eaf2a1..008f0f5da4 100644 --- a/src/source.f90 +++ b/src/source.f90 @@ -18,15 +18,17 @@ contains subroutine init_source() type(Particle), pointer :: p => null() - integer :: i,j,k - integer :: maxwork ! maxinum # of particles per processor - real(8) :: r(3) ! sampled coordinates - real(8) :: phi ! azimuthal angle - real(8) :: mu ! cosine of polar angle - real(8) :: E ! outgoing energy - real(8) :: p_min(3) ! minimum coordinates of source - real(8) :: p_max(3) ! maximum coordinates of source - character(250) :: msg ! error message + integer :: i ! loop index over processors + integer(8) :: j ! loop index over bank sites + integer :: k ! dummy loop index + integer :: maxwork ! maxinum # of particles per processor + real(8) :: r(3) ! sampled coordinates + real(8) :: phi ! azimuthal angle + real(8) :: mu ! cosine of polar angle + real(8) :: E ! outgoing energy + real(8) :: p_min(3) ! minimum coordinates of source + real(8) :: p_max(3) ! maximum coordinates of source + character(250) :: msg ! error message msg = 'Initializing source particles...' call message(msg, 6) diff --git a/src/tallies.f90 b/src/tallies.f90 index 3b728d6fbb..be3d6ff497 100644 --- a/src/tallies.f90 +++ b/src/tallies.f90 @@ -19,7 +19,7 @@ contains integer, intent(in) :: i_cycle ! index of current cycle - integer :: total_bank + integer(8) :: total_bank integer :: n integer :: ierr real(8) :: kcoll ! keff collision estimator @@ -40,7 +40,7 @@ contains #ifdef MPI ! Collect number bank sites onto master process - call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER, MPI_SUM, 0, & + call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & & MPI_COMM_WORLD, ierr) #else total_bank = n_bank diff --git a/src/types.f90 b/src/types.f90 index 48ca20af8f..ebfcb4c730 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -74,23 +74,23 @@ module types !===================================================================== type Particle - integer :: uid ! Unique ID - integer :: type ! Particle type (n, p, e, etc) - real(8) :: xyz(3) ! location - real(8) :: xyz_local(3) ! local location (after transformations) - real(8) :: uvw(3) ! directional cosines - real(8) :: wgt ! particle weight - real(8) :: E ! energy - integer :: IE ! index on energy grid - real(8) :: interp ! interpolation factor for energy grid - integer :: cell ! current cell - integer :: universe ! current universe - integer :: lattice ! current lattice - integer :: surface ! current surface - integer :: index_x ! lattice index for x direction - integer :: index_y ! lattice index for y direction - logical :: alive ! is particle alive? - integer :: n_coll ! # of collisions + integer(8) :: uid ! Unique ID + integer :: type ! Particle type (n, p, e, etc) + real(8) :: xyz(3) ! location + real(8) :: xyz_local(3) ! local location (after transformations) + real(8) :: uvw(3) ! directional cosines + real(8) :: wgt ! particle weight + real(8) :: E ! energy + integer :: IE ! index on energy grid + real(8) :: interp ! interpolation factor for energy grid + integer :: cell ! current cell + integer :: universe ! current universe + integer :: lattice ! current lattice + integer :: surface ! current surface + integer :: index_x ! lattice index for x direction + integer :: index_y ! lattice index for y direction + logical :: alive ! is particle alive? + integer :: n_coll ! # of collisions end type Particle !===================================================================== @@ -100,10 +100,10 @@ module types !===================================================================== type Bank - integer :: uid ! Unique ID - real(8) :: xyz(3) ! location of bank particle - real(8) :: uvw(3) ! diretional cosines - real(8) :: E ! energy + integer(8) :: uid ! Unique ID + real(8) :: xyz(3) ! location of bank particle + real(8) :: uvw(3) ! diretional cosines + real(8) :: E ! energy end type Bank !=====================================================================