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

(n,xn), inelastic scattering, Law 61. Parallel fission bank. Several bugs fixed.

Browse source
This commit is contained in:
Paul Romano 2011-03-21 14:52:49 +00:00
parent 6222890a8e
commit 095ce52afa
10 changed files with 697 additions and 92 deletions
Download patch file Download diff file Expand all files Collapse all files

52
ChangeLog
View file

@ -1,3 +1,55 @@
2011-03-21 Paul Romano <romano7@mit.edu>
* 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 <romano7@mit.edu>
* 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 <romano7@mit.edu>
* types.f90: The major changes are that the n_items and

4
src/Makefile
View file

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

13
src/ace.f90
View file

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

2
src/endf.f90
View file

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

1
src/global.f90
View file

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

8
src/input_sample
View file

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

21
src/main.f90
View file

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

150
src/mpi_routines.f90
View file

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

418
src/physics.f90
View file

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

120
src/source.f90
View file

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