crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

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.

Browse source
This commit is contained in:
Paul Romano 2011-04-21 20:24:08 +00:00
parent 5e18b4968c
commit 9a60b43946
11 changed files with 177 additions and 99 deletions
Download patch file Download diff file Expand all files Collapse all files

18
ChangeLog
View file

@ -1,3 +1,21 @@
2011-04-21 Paul Romano <romano7@mit.edu>
* 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 <romano7@mit.edu>
* main.f90: Simple timing for total time elapsed and bank

6
src/geometry.f90
View file

@ -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)

45
src/global.f90
View file

@ -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.

42
src/lattice_sample
View file

@ -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

12
src/main.f90
View file

@ -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

62
src/mpi_routines.f90
View file

@ -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)

6
src/output.f90
View file

@ -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)

19
src/physics.f90
View file

@ -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)

20
src/source.f90
View file

@ -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)

4
src/tallies.f90
View file

@ -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

42
src/types.f90
View file

@ -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
!=====================================================================