From 095ce52afad6d24c5ac1d97b81d4f00d1daebb68 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 21 Mar 2011 14:52:49 +0000 Subject: [PATCH] (n,xn), inelastic scattering, Law 61. Parallel fission bank. Several bugs fixed. --- ChangeLog | 52 ++++++ src/Makefile | 4 +- src/ace.f90 | 13 +- src/endf.f90 | 2 +- src/global.f90 | 1 + src/input_sample | 8 +- src/main.f90 | 21 ++- src/mpi_routines.f90 | 150 +++++++++++++++- src/physics.f90 | 418 ++++++++++++++++++++++++++++++++++++++----- src/source.f90 | 120 +++++++++---- 10 files changed, 697 insertions(+), 92 deletions(-) diff --git a/ChangeLog b/ChangeLog index 744001623..47de20a9c 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,55 @@ +2011-03-21 Paul Romano + + * ace.f90: Adjusted locators energy and angle distributions for + Law 61. + * global.f90: Added 'work' variable. + * main.f90: History loop is generalized. Added call to + synchronize_bank. + * mpi_routines.f90: Added fission bank synchronization from prior + code. Introduced MPI_bank type in setup_mpi. + * physics.f90: Added ability to simulate (n,xn), inelastic + scattering reactions. Further extended energy sampling to include + Law 61 (correlated energy-angle distributions). Fixed cases where + nu values outside of array for binary_search are treated + correctly. Fixed bugs in energy sampling (changed %energy to + %data). + * source.f90: Changed init_source to work in parallel. Minor + changes to get_source_particle. Added copy_from_bank for fission + bank synchronization. + +2011-03-18 Paul Romano + + * ace.f90: Added ability to read nu data through the read_nu_data + subroutine. Changed names of subroutines with 'parse' to + 'read'. Change 'size' variable to 'length' to prevent conflict + with SIZE intrinsic. Adjusted angle and energy distributions to + fit new types. Adjusted locator arrays within energy and angle + distribution data array. + * fileio.f90: Can now read verbosity from input file rather than + setting it arbitrarily in my source. + * geometry.f90: Set current pointers particle attributes in + find_cell. + * global.f90: Added keff, bank indices, parallel processing + variables, nu types. Moved default verbosity to set_defults. + * mpi_routines.f90: First addition on MPI parallelization + * output.f90: Added print_particle and print_reaction routines, + moved ou and eu to be module variables so that every routine + doesn't have to define them. + * physics.f90: Added n_fission and level_scatter + subroutines. Fixed bugs in sample_angle subroutine. Changed + rotate_angle to operate on particle. Added sample_energy + subroutine -- currently includes tabular equiprobable bins, + inelastic level scattering, continuous tabular distribution, + Maxwell spectrum, evaporation spectrum, Watt spectrum, and + Kalbach-Mann correlated scattering. + * source.f90: Added default values for source particles. Temporary + add_bank_sites routine written -- not sure if it will be used here + or elsewhere. + * types.f90: Added n_coll attribite to particle. Extended Bank + type to include direction and energy. Added AceDistAngle and + AceDistEnergy types and changed AceReaction accordingly. Changed + how nu data is stored in AceContinuous. + 2011-03-10 Paul Romano * types.f90: The major changes are that the n_items and diff --git a/src/Makefile b/src/Makefile index 84b0ea254..cbf9558a1 100644 --- a/src/Makefile +++ b/src/Makefile @@ -64,12 +64,12 @@ global.o: types.o main.o: global.o fileio.o output.o geometry.o mcnp_random.o \ source.o physics.o cross_section.o data_structures.o \ ace.o energy_grid.o -mpi_routines.o: global.o output.o types.o +mpi_routines.o: global.o output.o types.o mcnp_random.o source.o output.o: global.o types.o data_structures.o endf.o physics.o: types.o global.o mcnp_random.o geometry.o output.o \ search.o endf.o search.o: output.o -source.o: global.o mcnp_random.o +source.o: global.o mcnp_random.o output.o string.o: global.o output.o types.o: unittest.o: global.o energy_grid.o diff --git a/src/ace.f90 b/src/ace.f90 index fe858e22d..164b08c3b 100644 --- a/src/ace.f90 +++ b/src/ace.f90 @@ -584,7 +584,7 @@ contains integer :: NRa, NEa, NRb, NEb integer :: IDAT integer :: start, length, length_interp_data - integer :: i, j, k + integer :: i, j, k, l LED = JXS(10) LDIS = JXS(11) @@ -713,12 +713,23 @@ contains do j = 1,NE ! outgoing energy distribution NP = XSS(start + length + 2) + + ! adjust locators for angular distribution + do k = 1, NP + l = start + length + 2 + 3*NP + k + if (XSS(l) /= 0) XSS(l) = XSS(l) - LOCC - length_interp_data + end do + length = length + 2 + 4*NP do k = 1, NP ! outgoing angle distribution NP2 = XSS(start + length + 2) length = length + 2 + 3*NP2 end do + + ! adjust locators for energy distribution + k = start + 2 + 2*NR + NE + j + XSS(k) = XSS(k) - LOCC - length_interp_data end do case (66) diff --git a/src/endf.f90 b/src/endf.f90 index a76fea201..f232a39e2 100644 --- a/src/endf.f90 +++ b/src/endf.f90 @@ -51,7 +51,7 @@ contains case (33) string = '(n,nt)' case (34) - string = '(n,n He-3)' + string = '(n,nHe-3)' case (35) string = '(n,nd3a)' case (36) diff --git a/src/global.f90 b/src/global.f90 index 21332ddab..de856ed8e 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -59,6 +59,7 @@ module global 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 ! cycle keff real(8) :: keff diff --git a/src/input_sample b/src/input_sample index a90cdaa0b..df77b711f 100644 --- a/src/input_sample +++ b/src/input_sample @@ -3,14 +3,14 @@ cell 1 0 fill 37 -2 cell 100 37 40 -1 cell 101 37 41 1 -cell 2 0 40 2 -3 +cell 2 0 41 2 -3 surface 1 cylz 0 0 3 surface 2 cylz 0 0 5 -surface 3 cylz 0 0 8 +surface 3 cylz 0 0 7 material 40 -4.5 & - 92238.03c 1.0 + 92235.03c 1.0 material 41 -1.0 & 1001.03c 2.0 & @@ -20,5 +20,5 @@ source box -4 -4 -4 4 4 4 xs_library endfb7 xs_data /opt/serpent/xsdata/endfb7 -criticality 1 1 100 +criticality 1 1 10000 verbosity 8 diff --git a/src/main.f90 b/src/main.f90 index 0e7855b5c..287d98ffa 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -12,7 +12,8 @@ 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 + use mpi_routines, only: setup_mpi, synchronize_bank + use mpi implicit none @@ -95,16 +96,22 @@ contains subroutine run_problem() integer :: i, j - integer :: i_cycle - integer :: i_particle + integer :: ierr + integer :: i_cycle ! cycle index + integer :: i_particle ! history index + integer :: total_bank ! total number of particles banked type(Particle), pointer :: p => null() + character(250) :: msg + + msg = "Running problem..." + call message(msg, 6) CYCLE_LOOP: do i_cycle = 1, n_cycles ! Set all tallies to zero n_bank = 0 - HISTORY_LOOP: do j = 1, n_particles + HISTORY_LOOP: do ! grab source particle from bank p => get_source_particle() @@ -114,7 +121,7 @@ contains end if ! set random number seed - i_particle = (i-1)*n_particles + j + i_particle = (i_cycle-1)*n_particles + p % uid call RN_init_particle(int(i_particle,8)) ! transport particle @@ -122,8 +129,10 @@ contains end do HISTORY_LOOP + call RN_init_particle(int(i_cycle,8)) + call synchronize_bank() + ! Collect results and statistics - print *, n_bank ! print cycle information diff --git a/src/mpi_routines.f90 b/src/mpi_routines.f90 index 7eae9e6e3..350ba0e86 100644 --- a/src/mpi_routines.f90 +++ b/src/mpi_routines.f90 @@ -3,6 +3,8 @@ module mpi_routines use mpi use global use output, only: message, error + use mcnp_random, only: rang + use source, only: copy_from_bank implicit none @@ -19,8 +21,11 @@ contains subroutine setup_mpi() - integer :: ierr - character(250) :: msg + integer :: ierr ! Error status + integer :: bank_blocks(2) ! Count for each datatype + integer :: bank_types(2) ! Datatypes + integer(kind=8) :: bank_disp(2) ! Displacements + character(250) :: msg ! Error message ! Initialize MPI call MPI_INIT(ierr) @@ -51,6 +56,147 @@ contains master = .false. end if + ! Define MPI_bank for fission sites + bank_blocks = (/ 1, 7 /) + bank_disp = (/ 0_8, 4_8 /) + bank_types = (/ MPI_INTEGER, MPI_REAL8 /) + call MPI_TYPE_CREATE_STRUCT(2, bank_blocks, bank_disp, & + & bank_types, MPI_bank, ierr) + call MPI_TYPE_COMMIT(MPI_bank, ierr) + end subroutine setup_mpi + +!===================================================================== +! SYNCHRONIZE_BANK +!===================================================================== + + subroutine synchronize_bank() + + implicit none + + integer :: i, j, k ! loop indices + integer :: status(MPI_STATUS_SIZE) ! message status + 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 + 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 + + ! 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, & + & MPI_COMM_WORLD, ierr) + finish = start + n_bank + total = finish + call MPI_BCAST(total, 1, MPI_INTEGER, n_procs - 1, & + & MPI_COMM_WORLD, ierr) + + allocate(temp_sites(2*work)) + count = 0 ! Index for source_bank + index = 0 ! Index for uid -- must account for all nodes + + if (total < n_particles) then + sites_needed = mod(n_particles,total) + else + sites_needed = n_particles + end if + + ! ================================================================ + ! SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES + do i = 1, total + + ! If there are less than n_particles particles banked, automatically + ! add int(n_particles/total) sites to temp_sites + if (total < n_particles) then + do j = 1,int(n_particles/total) + index = index + 1 + ! If index is within this node's range, add site to source + if (i > start .and. i <= finish) then + count = count + 1 + temp_sites(count) = fission_bank(i-start) + temp_sites(count) % uid = index + end if + end do + end if + + ! Randomly sample sites needed + sites_remaining = total - i + 1 + if (sites_needed == sites_remaining .or. & + rang() < real(sites_needed)/real(sites_remaining)) then + index = index + 1 + if (i > start .and. i <= finish) then + count = count + 1 + temp_sites(count) = fission_bank(i-start) + temp_sites(count) % uid = index + end if + sites_needed = sites_needed-1 + end if + end do + + ! Determine how many sites to send to adjacent nodes + send_to_left = bank_first - temp_sites(1)%uid + send_to_right = temp_sites(1)%uid + (count-1) - bank_last + + ! ================================================================ + ! 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 = count - send_to_right + 1 + call MPI_SEND(temp_sites(i), send_to_right, MPI_bank, rank+1, 0, & + & MPI_COMM_WORLD, ierr) + else if (send_to_right < 0) then + call MPI_RECV(right_bank, -send_to_right, MPI_bank, rank+1, 1, & + & MPI_COMM_WORLD, status, ierr) + end if + + if (send_to_left < 0) then + call MPI_RECV(left_bank, -send_to_left, MPI_bank, rank-1, 0, & + & MPI_COMM_WORLD, status, ierr) + else if (send_to_left > 0) then + call MPI_SEND(temp_sites(1), send_to_left, MPI_bank, rank-1, 1, & + & MPI_COMM_WORLD, ierr) + end if + + ! ================================================================ + ! RECONSTRUCT SOURCE BANK + if (send_to_left < 0 .and. send_to_right >= 0) then + i = -send_to_left ! size of first block + j = count - send_to_right ! size of second block + call copy_from_bank(left_bank, 1, i) + call copy_from_bank(temp_sites, i+1, j) + else if (send_to_left >= 0 .and. send_to_right < 0) then + i = count - 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) + call copy_from_bank(right_bank, i+1, j) + else if (send_to_left >= 0 .and. send_to_right >= 0) then + i = count - 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 = count + k = -send_to_right + call copy_from_bank(left_bank, 1, i) + call copy_from_bank(temp_sites, i+1, j) + call copy_from_bank(right_bank, i+j+1, k) + end if + + deallocate(left_bank ) + deallocate(right_bank) + deallocate(temp_sites) + + end subroutine synchronize_bank end module mpi_routines diff --git a/src/physics.f90 b/src/physics.f90 index 0d9c3edd6..47aa44a1f 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -4,7 +4,7 @@ module physics use geometry, only: find_cell, dist_to_boundary, cross_boundary use types, only: Particle use mcnp_random, only: rang - use output, only: error, message + use output, only: error, message, print_particle use search, only: binary_search use endf, only: reaction_name @@ -213,6 +213,10 @@ contains select case (rxn%MT) case (2) call elastic_scatter(p, table%awr) + case (11, 16, 17, 24, 25, 30, 37, 41, 42) + call n_xn(p, table, rxn) + case (22, 23, 28, 29, 32:36, 43, 44, 91) + call inelastic_scatter(p, table, rxn) case (18:21, 38) call n_fission(p, table, rxn) case (51:90) @@ -359,9 +363,17 @@ contains ! do binary search over tabuled energies to determine ! appropriate index and interpolation factor - j = binary_search(table % nu_t_data(loc+1), NE, E) - f = (E - table % nu_t_data(loc+j)) / & - & (table % nu_t_data(loc+j+1) - table % nu_t_data(loc+j)) + if (E < table % nu_t_data(loc+1)) then + j = 1 + f = 0.0 + elseif (E > table % nu_t_data(loc+NE)) then + j = NE - 1 + f = 1.0 + else + j = binary_search(table % nu_t_data(loc+1), NE, E) + f = (E - table % nu_t_data(loc+j)) / & + & (table % nu_t_data(loc+j+1) - table % nu_t_data(loc+j)) + end if ! determine nu total loc = loc + NE @@ -400,9 +412,17 @@ contains ! do binary search over tabuled energies to determine ! appropriate index and interpolation factor - j = binary_search(table % nu_p_data(loc+1), NE, E) - f = (E - table % nu_p_data(loc+j)) / & - & (table % nu_p_data(loc+j+1) - table % nu_p_data(loc+j)) + if (E < table % nu_p_data(loc+1)) then + j = 1 + f = 0.0 + elseif (E > table % nu_p_data(loc+NE)) then + j = NE - 1 + f = 1.0 + else + j = binary_search(table % nu_p_data(loc+1), NE, E) + f = (E - table % nu_p_data(loc+j)) / & + & (table % nu_p_data(loc+j+1) - table % nu_p_data(loc+j)) + end if ! determine nu total loc = loc + NE @@ -429,9 +449,17 @@ contains ! do binary search over tabuled energies to determine ! appropriate index and interpolation factor - j = binary_search(table % nu_d_data(loc+1), NE, E) - f = (E - table % nu_d_data(loc+j)) / & - & (table % nu_d_data(loc+j+1) - table % nu_d_data(loc+j)) + if (E < table % nu_d_data(loc+1)) then + j = 1 + f = 0.0 + elseif (E > table % nu_d_data(loc+NE)) then + j = NE - 1 + f = 1.0 + else + j = binary_search(table % nu_d_data(loc+1), NE, E) + f = (E - table % nu_d_data(loc+j)) / & + & (table % nu_d_data(loc+j+1) - table % nu_d_data(loc+j)) + end if ! determine nu total loc = loc + NE @@ -540,6 +568,140 @@ contains end subroutine level_scatter +!===================================================================== +! INELASTIC_SCATTER +!===================================================================== + + subroutine inelastic_scatter(p, table, rxn) + + type(Particle), pointer :: p + type(AceContinuous), pointer :: table + type(AceReaction), pointer :: rxn + + integer :: law + real(8) :: A ! atomic weight ratio of nuclide + real(8) :: E_in ! incoming energy + real(8) :: mu ! cosine of scattering angle + real(8) :: E ! outgoing energy in laboratory + real(8) :: E_cm ! outgoing energy in center-of-mass + character(250) :: msg ! error message + + ! copy energy of neutron + E_in = p % E + + ! determine A + A = table % awr + + ! determine secondary energy distribution law + law = rxn % edist % law + + ! sample scattering angle + if (rxn % has_angle_dist) then + mu = sample_angle(rxn, E_in) + end if + + ! sample outgoing energy + if (law == 44 .or. law == 61) then + call sample_energy(rxn, E_in, E, mu) + else + call sample_energy(rxn, E_in, E) + end if + + ! if scattering system is in center-of-mass, transfer cosine of + ! scattering angle and outgoing energy from CM to LAB + if (rxn % TY < 0) then + E_cm = E + + ! determine outgoing energy in lab + E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) & + & / ((A+ONE)*(A+ONE)) + + ! determine outgoing angle in lab + mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E) + end if + + ! change direction of particle + call rotate_angle(p, mu) + + ! change energy of particle + p % E = E + + ! find energy index, interpolation factor + call find_energy_index(p) + + end subroutine inelastic_scatter + +!===================================================================== +! N_XN +!===================================================================== + + subroutine n_xn(p, table, rxn) + + type(Particle), pointer :: p + type(AceContinuous), pointer :: table + type(AceReaction), pointer :: rxn + + integer :: i ! loop index + integer :: n_secondary ! number of secondary particles + integer :: law + real(8) :: A ! atomic weight ratio of nuclide + real(8) :: E_in ! incoming energy + real(8) :: mu ! cosine of scattering angle + real(8) :: E ! outgoing energy in laboratory + real(8) :: E_cm ! outgoing energy in center-of-mass + character(250) :: msg ! error message + + ! copy energy of neutron + E_in = p % E + + ! determine A + A = table % awr + + ! determine secondary energy distribution law + law = rxn % edist % law + + ! sample scattering angle + if (rxn % has_angle_dist) then + mu = sample_angle(rxn, E_in) + end if + + ! sample outgoing energy + if (law == 44 .or. law == 61) then + call sample_energy(rxn, E_in, E, mu) + else + call sample_energy(rxn, E_in, E) + end if + + ! if scattering system is in center-of-mass, transfer cosine of + ! scattering angle and outgoing energy from CM to LAB + if (rxn % TY < 0) then + E_cm = E + + ! determine outgoing energy in lab + E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) & + & / ((A+ONE)*(A+ONE)) + + ! determine outgoing angle in lab + mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E) + end if + + ! change direction of particle + call rotate_angle(p, mu) + + ! change energy of particle + p % E = E + + ! find energy index, interpolation factor + call find_energy_index(p) + + ! produce secondary particles in source bank + n_secondary = abs(rxn % TY) - 1 + do i = 1, n_secondary + ! create particle in source bank + end do + + end subroutine n_xn + !===================================================================== ! N_GAMMA !===================================================================== @@ -735,6 +897,7 @@ contains integer :: NET ! number of outgoing energies integer :: INTTp ! combination of INTT and ND integer :: INTT ! 1 = histogram, 2 = linear-linear + integer :: JJ ! 1 = histogram, 2 = linear-linear integer :: ND ! number of discrete lines integer :: NP ! number of points in distribution @@ -753,6 +916,11 @@ contains real(8) :: Watt_a, Watt_b ! Watt spectrum parameters + real(8) :: mu_k ! angular cosine in bin k + real(8) :: mu_k1 ! angular cosine in bin k+1 + real(8) :: p_k ! angular pdf in bin k + real(8) :: p_k1 ! angular pdf in bin k+1 + real(8) :: E_cm real(8) :: xi1, xi2, xi3, xi4 real(8) :: r ! interpolation factor on incoming energy @@ -843,13 +1011,13 @@ contains if (E_in < rxn % edist % data(loc+1)) then i = 1 r = 0.0 - elseif (E_in > rxn % edist % energy(loc+NE)) then + elseif (E_in > rxn % edist % data(loc+NE)) then i = NE - 1 r = 1.0 else - i = binary_search(rxn % edist % energy(loc+1), NE, E_in) - r = (E_in - rxn%edist%energy(loc+i)) / & - & (rxn%edist%energy(loc+i+1) - rxn%edist%energy(loc+i)) + i = binary_search(rxn % edist % data(loc+1), NE, E_in) + r = (E_in - rxn%edist%data(loc+i)) / & + & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) end if ! Sample between the ith and (i+1)th bin @@ -861,11 +1029,11 @@ contains end if ! interpolation on grid i for energy E1 - loc = rxn % edist % data(2 + 2*NR + NE + i) + 2 ! start of EOUT(i) + loc = rxn%edist%data(2 + 2*NR + NE + i) + 2 ! start of EOUT(i) E_i_1 = rxn%edist%data(loc + 1) E_i_K = rxn%edist%data(loc + NE) - loc = rxn % edist % data(2 + 2*NR + NE + i + 1) + 2 ! start of EOUT(i+1) + loc = rxn%edist%data(2 + 2*NR + NE + i + 1) + 2 ! start of EOUT(i+1) E_i1_1 = rxn%edist%data(loc + 1) E_i1_K = rxn%edist%data(loc + NE) @@ -927,7 +1095,7 @@ contains call error(msg) end if - ! Now interpolate between incident enregy bins i and i + 1 + ! Now interpolate between incident energy bins i and i + 1 if (l == i) then E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) else @@ -956,13 +1124,13 @@ contains if (E_in < rxn % edist % data(loc+1)) then i = 1 r = 0.0 - elseif (E_in > rxn % edist % energy(loc+NE)) then + elseif (E_in > rxn % edist % data(loc+NE)) then i = NE - 1 r = 1.0 else - i = binary_search(rxn % edist % energy(loc+1), NE, E_in) - r = (E_in - rxn%edist%energy(loc+i)) / & - & (rxn%edist%energy(loc+i+1) - rxn%edist%energy(loc+i)) + i = binary_search(rxn % edist % data(loc+1), NE, E_in) + r = (E_in - rxn%edist%data(loc+i)) / & + & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) end if ! determine nuclear temperature from tabulated function @@ -993,13 +1161,13 @@ contains if (E_in < rxn % edist % data(loc+1)) then i = 1 r = 0.0 - elseif (E_in > rxn % edist % energy(loc+NE)) then + elseif (E_in > rxn % edist % data(loc+NE)) then i = NE - 1 r = 1.0 else - i = binary_search(rxn % edist % energy(loc+1), NE, E_in) - r = (E_in - rxn%edist%energy(loc+i)) / & - & (rxn%edist%energy(loc+i+1) - rxn%edist%energy(loc+i)) + i = binary_search(rxn % edist % data(loc+1), NE, E_in) + r = (E_in - rxn%edist%data(loc+i)) / & + & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) end if ! determine nuclear temperature from tabulated function @@ -1035,13 +1203,13 @@ contains if (E_in < rxn % edist % data(loc+1)) then i = 1 r = 0.0 - elseif (E_in > rxn % edist % energy(loc+NE)) then + elseif (E_in > rxn % edist % data(loc+NE)) then i = NE - 1 r = 1.0 else - i = binary_search(rxn % edist % energy(loc+1), NE, E_in) - r = (E_in - rxn%edist%energy(loc+i)) / & - & (rxn%edist%energy(loc+i+1) - rxn%edist%energy(loc+i)) + i = binary_search(rxn % edist % data(loc+1), NE, E_in) + r = (E_in - rxn%edist%data(loc+i)) / & + & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) end if ! determine Watt parameter 'a' from tabulated function @@ -1065,13 +1233,13 @@ contains if (E_in < rxn % edist % data(loc+1)) then i = 1 r = 0.0 - elseif (E_in > rxn % edist % energy(loc+NE)) then + elseif (E_in > rxn % edist % data(loc+NE)) then i = NE - 1 r = 1.0 else - i = binary_search(rxn % edist % energy(loc+1), NE, E_in) - r = (E_in - rxn%edist%energy(loc+i)) / & - & (rxn%edist%energy(loc+i+1) - rxn%edist%energy(loc+i)) + i = binary_search(rxn % edist % data(loc+1), NE, E_in) + r = (E_in - rxn%edist%data(loc+i)) / & + & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) end if ! determine Watt parameter 'b' from tabulated function @@ -1107,13 +1275,13 @@ contains if (E_in < rxn % edist % data(loc+1)) then i = 1 r = 0.0 - elseif (E_in > rxn % edist % energy(loc+NE)) then + elseif (E_in > rxn % edist % data(loc+NE)) then i = NE - 1 r = 1.0 else - i = binary_search(rxn % edist % energy(loc+1), NE, E_in) - r = (E_in - rxn%edist%energy(loc+i)) / & - & (rxn%edist%energy(loc+i+1) - rxn%edist%energy(loc+i)) + i = binary_search(rxn % edist % data(loc+1), NE, E_in) + r = (E_in - rxn%edist%data(loc+i)) / & + & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) end if ! Sample between the ith and (i+1)th bin @@ -1124,14 +1292,19 @@ contains l = i end if - ! interpolation on grid i for energy E1 - loc = rxn % edist % data(2 + 2*NR + NE + i) + 2 ! start of EOUT(i) + ! determine endpoints on grid i + loc = rxn%edist%data(2+2*NR+NE + i) ! start of LDAT for i + NP = rxn%edist%data(loc + 2) + loc = loc + 2 E_i_1 = rxn%edist%data(loc + 1) - E_i_K = rxn%edist%data(loc + NE) + E_i_K = rxn%edist%data(loc + NP) - loc = rxn % edist % data(2 + 2*NR + NE + i + 1) + 2 ! start of EOUT(i+1) + ! determine endpoints on grid i+1 + loc = rxn%edist%data(2+2*NR+NE + i+1) ! start of LDAT for i+1 + NP = rxn%edist%data(loc + 2) + loc = loc + 2 E_i1_1 = rxn%edist%data(loc + 1) - E_i1_K = rxn%edist%data(loc + NE) + E_i1_K = rxn%edist%data(loc + NP) E_1 = E_i_1 + r*(E_i1_1 - E_i_1) E_K = E_i_K + r*(E_i1_K - E_i_K) @@ -1205,7 +1378,7 @@ contains call error(msg) end if - ! Now interpolate between incident enregy bins i and i + 1 + ! Now interpolate between incident energy bins i and i + 1 if (l == i) then E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) else @@ -1231,6 +1404,165 @@ contains call error(msg) end if + ! read number of interpolation regions and incoming energies + NR = rxn % edist % data(1) + NE = rxn % edist % data(2 + 2*NR) + if (NR > 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to sample correlated energy-angle distribution." + call error(msg) + end if + + ! find energy bin and calculate interpolation factor -- if the + ! energy is outside the range of the tabulated energies, choose + ! the first or last bins + loc = 2 + 2*NR + if (E_in < rxn % edist % data(loc+1)) then + i = 1 + r = 0.0 + elseif (E_in > rxn % edist % data(loc+NE)) then + i = NE - 1 + r = 1.0 + else + i = binary_search(rxn % edist % data(loc+1), NE, E_in) + r = (E_in - rxn%edist%data(loc+i)) / & + & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + end if + + ! Sample between the ith and (i+1)th bin + xi2 = rang() + if (r > xi2) then + l = i + 1 + else + l = i + end if + + ! determine endpoints on grid i + loc = rxn%edist%data(2+2*NR+NE + i) ! start of LDAT for i + NP = rxn%edist%data(loc + 2) + loc = loc + 2 + E_i_1 = rxn%edist%data(loc + 1) + E_i_K = rxn%edist%data(loc + NP) + + ! determine endpoints on grid i+1 + loc = rxn%edist%data(2+2*NR+NE + i+1) ! start of LDAT for i+1 + NP = rxn%edist%data(loc + 2) + loc = loc + 2 + E_i1_1 = rxn%edist%data(loc + 1) + E_i1_K = rxn%edist%data(loc + NP) + + E_1 = E_i_1 + r*(E_i1_1 - E_i_1) + E_K = E_i_K + r*(E_i1_K - E_i_K) + + ! determine location of outgoing energies, pdf, cdf for E(l) + loc = rxn % edist % data(2 + 2*NR + NE + l) + + ! determine type of interpolation and number of discrete lines + INTTp = rxn % edist % data(loc + 1) + NP = rxn % edist % data(loc + 2) + if (INTTp > 10) then + INTT = mod(INTTp,10) + ND = (INTTp - INTT)/10 + else + INTT = INTTp + ND = 0 + end if + + if (ND > 0) then + ! discrete lines present + msg = "Discrete lines in continuous tabular distributed not & + &yet supported" + call error(msg) + end if + + ! determine outgoing energy bin + xi1 = rang() + loc = loc + 2 ! start of EOUT + c_k = rxn % edist % data(loc + 2*NP + 1) + do k = 1, NP-1 + c_k1 = rxn % edist % data(loc + 2*NP + k+1) + if (xi1 < c_k1) exit + c_k = c_k1 + end do + + E_l_k = rxn % edist % data(loc+k) + p_l_k = rxn % edist % data(loc+NP+k) + if (INTT == HISTOGRAM) then + ! Histogram interpolation + E_out = E_l_k + (xi1 - c_k)/p_l_k + + elseif (INTT == LINEARLINEAR) then + ! Linear-linear interpolation -- not sure how you come about + ! the formula given in the MCNP manual + E_l_k1 = rxn % edist % data(loc+k+1) + p_l_k1 = rxn % edist % data(loc+NP+k+1) + + ! Find E prime + frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) + if (frac == ZERO) then + E_out = E_l_k + (xi1 - c_k)/p_l_k + else + E_out = E_l_k + (sqrt(p_l_k*p_l_k + 2*frac*(xi1 - c_k)) - & + & p_l_k)/frac + end if + else + msg = "Unknown interpolation type: " // trim(int_to_str(INTT)) + call error(msg) + end if + + ! Now interpolate between incident energy bins i and i + 1 + if (l == i) then + E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) + else + E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) + end if + + ! Find location of correlated angular distribution + loc = rxn % edist % data(loc+3*NP+k) + + ! Check if angular distribution is isotropic + if (loc == 0) then + mu_out = TWO * rang() - ONE + return + end if + + ! interpolation type and number of points in angular distribution + JJ = rxn % edist % data(loc + 1) + NP = rxn % edist % data(loc + 2) + + ! determine outgoing cosine bin + xi3 = rang() + loc = loc + 2 + c_k = rxn % edist % data(loc + 2*NP + 1) + do k = 1, NP-1 + c_k1 = rxn % edist % data(loc + 2*NP + k+1) + if (xi3 < c_k1) exit + c_k = c_k1 + end do + + p_k = rxn % edist % data(loc + NP + k) + mu_k = rxn % edist % data(loc + k) + if (JJ == HISTOGRAM) then + ! Histogram interpolation + mu_out = mu_k + (xi3 - c_k)/p_k + + elseif (JJ == LINEARLINEAR) then + ! Linear-linear interpolation -- not sure how you come about + ! the formula given in the MCNP manual + p_k1 = rxn % edist % data(loc + NP + k+1) + mu_k1 = rxn % edist % data(loc + k+1) + + frac = (p_k1 - p_k)/(mu_k1 - mu_k) + if (frac == ZERO) then + mu_out = mu_k + (xi3 - c_k)/p_k + else + mu_out = mu_k + (sqrt(p_k*p_k + 2*frac*(xi3 - c_k))-p_k)/frac + end if + else + msg = "Unknown interpolation type: " // trim(int_to_str(JJ)) + call error(msg) + end if + case (66) ! ============================================================= ! N-BODY PHASE SPACE DISTRIBUTION diff --git a/src/source.f90 b/src/source.f90 index 15c6c3386..e1d83208d 100644 --- a/src/source.f90 +++ b/src/source.f90 @@ -2,6 +2,7 @@ module source use global use mcnp_random, only: rang, RN_init_particle + use output, only: message implicit none @@ -15,48 +16,70 @@ contains subroutine init_source() - integer :: i,j - real(8) :: r(3) - real(8) :: phi ! azimuthal angle - real(8) :: mu ! cosine of polar angle - real(8) :: p_min(3), p_max(3) 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) :: 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) ! Allocate fission and source banks - allocate( source_bank(n_particles) ) - allocate( fission_bank(3*n_particles) ) + allocate(source_bank(n_particles)) + allocate(fission_bank(3*n_particles)) + + ! Determine maximum amount of particles to simulate on each + ! processor + maxwork = ceiling(real(n_particles)/n_procs) ! Check external source type - if ( external_source%type == SRC_BOX ) then + if (external_source%type == SRC_BOX) then p_min = external_source%values(1:3) p_max = external_source%values(4:6) end if ! Initialize first cycle source bank - do i = 1, n_particles - call RN_init_particle(int(i,8)) - p => source_bank(i) + do i = 0, n_procs - 1 + if (rank == i) then + ! UID's of first and last source particles + bank_first = i*maxwork + 1 + bank_last = min((i+1)*maxwork, n_particles) - ! position - r = (/ (rang(), j = 1,3) /) - p % uid = i - p % xyz = p_min + r*(p_max - p_min) + ! number of particles for this processor + work = bank_last - bank_first + 1 - ! angle - phi = TWO*PI*rang() - mu = TWO*rang() - ONE - p % uvw(1) = mu - p % uvw(2) = sqrt(ONE - mu*mu) * cos(phi) - p % uvw(3) = sqrt(ONE - mu*mu) * sin(phi) + do j = bank_first, bank_last + p => source_bank(j - bank_first + 1) - ! set defaults - p % type = NEUTRON - p % cell = 0 - p % surface = 0 - p % universe = 0 - p % wgt = ONE - p % alive = .true. + ! initialize random number seed + call RN_init_particle(int(j,8)) + ! position + r = (/ (rang(), k = 1,3) /) + p % uid = j + p % xyz = p_min + r*(p_max - p_min) + + ! angle + phi = TWO*PI*rang() + mu = TWO*rang() - ONE + p % uvw(1) = mu + p % uvw(2) = sqrt(ONE - mu*mu) * cos(phi) + p % uvw(3) = sqrt(ONE - mu*mu) * sin(phi) + + ! set defaults + p % type = NEUTRON + p % cell = 0 + p % surface = 0 + p % universe = 0 + p % wgt = ONE + p % alive = .true. + end do + end if end do ! Reset source index @@ -68,21 +91,21 @@ contains ! GET_SOURCE_PARTICLE returns the next source particle !===================================================================== - function get_source_particle() + function get_source_particle() result(p) - type(Particle), pointer :: get_source_particle + type(Particle), pointer :: p ! increment index source_index = source_index + 1 ! if at end of bank, return null pointer - if (source_index > size(source_bank)) then - get_source_particle => null() + if (source_index > work) then + p => null() return end if ! point to next source particle - get_source_particle => source_bank(source_index) + p => source_bank(source_index) end function get_source_particle @@ -110,4 +133,35 @@ contains end subroutine add_bank_sites +!===================================================================== +! COPY_FROM_BANK +!===================================================================== + + subroutine copy_from_bank(temp_bank, index, n_sites) + + type(Bank), intent(in) :: temp_bank(n_sites) + integer, intent(in) :: index ! starting index in source_bank + integer, intent(in) :: n_sites ! # of bank sites to copy + + integer :: i ! index in temp_bank + integer :: j ! index in source_bank + + do i = 1, n_sites + j = index + i - 1 + source_bank(j) % uid = temp_bank(i) % uid + source_bank(j) % xyz = temp_bank(i) % xyz + source_bank(j) % uvw = temp_bank(i) % uvw + source_bank(j) % E = temp_bank(i) % E + + ! set defaults + source_bank(j) % type = NEUTRON + source_bank(j) % cell = 0 + source_bank(j) % surface = 0 + source_bank(j) % universe = 0 + source_bank(j) % wgt = ONE + source_bank(j) % alive = .true. + end do + + end subroutine copy_from_bank + end module source