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

Fixed major bugs, namely -- material total xs, MPI_BANK, evaporation spectrum, interpolation. Added initial keff ability. Added Jezebel benchmark.

Browse source
This commit is contained in:
Paul Romano 2011-03-22 14:55:26 +00:00
parent 095ce52afa
commit eee77ec2f3
12 changed files with 328 additions and 169 deletions
Download patch file Download diff file Expand all files Collapse all files

22
ChangeLog
View file

@ -1,3 +1,25 @@
2011-03-22 Paul Romano <romano7@mit.edu>
* cross_section.f90: Calculation of material total cross section
was completely wrong -- fixed it up.
* global.f90: Added default keff.
* mpi_routines.f90: Fixed major bug in how MPI_BANK type was
created. Solution should be machine independent. Parallel fission
bank now works.
* physics.f90: Combined level_scatter into inelastic_scatter since
the physics is the same. Simplifying sampling of angle (isotropic
cases handled inside sample_angle). Fixed interpolation in many
places (f and (1-f) were swapped). Fixed distance to
collision. Fixed major bug in evaporation spectrum (restriction
energy was wrong).
* source.f90: Moved initial sampling of Watt spectrum for first
cycle particles to init_source inside of being in physics module.
* tallies.f90: Initial addition of tallies module. Only has one
routine, calculate_keff, that is used to do end-of-cycle keff
calculations. This will probably change significantly.
* unittest.f90: Added test_mpi_bank useful for debugging the
MPI_BANK type.
2011-03-21 Paul Romano <romano7@mit.edu>
* ace.f90: Adjusted locators energy and angle distributions for

24
benchmarks/pu-met-fast-001 Normal file
View file

@ -0,0 +1,24 @@
# ====================================================================
# Description: Bare Pu-239 Jezebel Benchmark
# Case: PU-MET-FAST-001
# Written By: Paul Romano
# Date: 3/21/2011
# ====================================================================
# cell no. universe material surfaces
cell 1 0 1 -1
surface 1 sph 0. 0. 0. 6.3849
source_box -4 -4 -4 4 4 4
material 1 0.04029014 &
94239.03c 0.037047 &
94240.03c 0.0017512 &
94241.03c 0.00011674 &
31000.03c 0.0013752
xs_library endfb7
xs_data /opt/serpent/xsdata/endfb7
criticality 150 25 1000
verbosity 7

10
src/Makefile
View file

@ -17,6 +17,7 @@ modules = ace.f90 \
search.f90 \
source.f90 \
string.f90 \
tallies.f90 \
types.f90
main_objects = $(modules:.f90=.o) $(main:.f90=.o)
@ -41,7 +42,7 @@ clean:
neat:
@rm -f *.o *.mod
unittest: $(test_objects)
$(F90) $(F90FLAGS) $(MPI_COMPILE_FLAGS) $(test_objects) -o $@
$(MPIF90) -f90=$(F90) $(MPI_LINK_FLAGS) $(F90FLAGS) $(test_objects) -o $@
#--------------------------------------------------------------------
# Rules & misc
@ -63,13 +64,14 @@ geometry.o: types.o global.o output.o string.o data_structures.o
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
ace.o energy_grid.o tallies.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 output.o
source.o: global.o mcnp_random.o output.o physics.o
string.o: global.o output.o
tallies.o: global.o output.o
types.o:
unittest.o: global.o energy_grid.o
unittest.o: global.o energy_grid.o mpi_routines.o

57
src/cross_section.f90
View file

@ -125,12 +125,20 @@ contains
subroutine material_total_xs()
type(Material), pointer :: mat => null()
type(AceContinuous), pointer :: acecont => null()
type(Material), pointer :: mat => null()
type(AceContinuous), pointer :: table => null()
integer :: i, j, k
integer :: index
integer :: IE
real(8) :: xs
real(8) :: density
real(8) :: density_i
real(8) :: val
real(8) :: r
real(8) :: E_i
real(8) :: E_i1
real(8) :: sigma_i
real(8) :: sigma_i1
character(250) :: msg
msg = "Creating material total cross-sections..."
@ -143,21 +151,40 @@ contains
mat => materials(i)
density = mat%atom_density
allocate(mat%total_xs(n_grid))
! loop over each energy on grid
do j = 1, n_grid
mat%total_xs(j) = 0.0_8
! loop over each isotope in material
do k = 1, mat%n_isotopes
! find ACE cross section table for this isotope
index = mat%table(k)
acecont => xs_continuous(index)
! initialize total cross-section
mat % total_xs = ZERO
! loop over each isotope
do j = 1, mat%n_isotopes
index = mat % table(j)
table => xs_continuous(index)
! find atom density of isotope
density_i = density * mat % atom_percent(j)
! loop over points in union energy grid
do k = 1, n_grid
! find nuclide grid index
IE = table % grid_index(k)
! get nuclide grid energy and cross-section value
E_i = table % energy(IE)
sigma_i = table % sigma_t(IE)
! interpolate as necessary
if (E_i < e_grid(k)) then
E_i1 = table % energy(IE + 1)
sigma_i1 = table % sigma_t(IE + 1)
r = (e_grid(k) - E_i)/(E_i1 - E_i)
val = (ONE - r)*sigma_i + r*sigma_i1
else
val = sigma_i
end if
! add contribution to total cross-section
mat%total_xs(j) = mat%total_xs(j) + density * &
mat%atom_percent(j) * acecont%sigma_t(index)
mat%total_xs(k) = mat%total_xs(k) + density_i * val
end do
end do
end do

3
src/global.f90
View file

@ -200,6 +200,9 @@ contains
! Default verbosity
verbosity = 5
! Defualt multiplication factor
keff = ONE
end subroutine set_defaults
!=====================================================================

2
src/input_sample
View file

@ -20,5 +20,5 @@ source box -4 -4 -4 4 4 4
xs_library endfb7
xs_data /opt/serpent/xsdata/endfb7
criticality 1 1 10000
criticality 15 5 10000
verbosity 8

10
src/main.f90
View file

@ -4,7 +4,7 @@ program main
use fileio, only: read_input, read_command_line, read_count, &
& normalize_ao, build_universe
use output, only: title, echo_input, message, warning, error, &
& print_summary
& print_summary, print_particle
use geometry, only: sense, cell_contains, neighbor_lists
use mcnp_random, only: RN_init_problem, rang, RN_init_particle
use source, only: init_source, get_source_particle
@ -13,6 +13,7 @@ program main
use ace, only: read_xs
use energy_grid, only: unionized_grid, original_indices
use mpi_routines, only: setup_mpi, synchronize_bank
use tallies, only: calculate_keff
use mpi
implicit none
@ -72,15 +73,15 @@ program main
call material_total_xs()
if (master) then
! call echo_input()
! call print_summary()
call echo_input()
call print_summary()
end if
! create source particles
call init_source()
! start problem
surfaces(3)%bc = BC_VACUUM
surfaces(1)%bc = BC_VACUUM
call run_problem()
! deallocate arrays
@ -133,6 +134,7 @@ contains
call synchronize_bank()
! Collect results and statistics
call calculate_keff(i_cycle)
! print cycle information

66
src/mpi_routines.f90
View file

@ -4,11 +4,11 @@ module mpi_routines
use global
use output, only: message, error
use mcnp_random, only: rang
use source, only: copy_from_bank
use source, only: copy_from_bank, source_index
implicit none
integer :: MPI_bank ! MPI datatype for fission bank
integer :: MPI_BANK ! MPI datatype for fission bank
integer(8) :: bank_index ! Fission bank site unique identifier
contains
@ -21,11 +21,14 @@ contains
subroutine setup_mpi()
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
integer :: i
integer :: ierr ! Error status
integer :: bank_blocks(4) ! Count for each datatype
integer :: bank_types(4) ! Datatypes
integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements
integer(MPI_ADDRESS_KIND) :: base
character(250) :: msg ! Error message
type(Bank) :: b
! Initialize MPI
call MPI_INIT(ierr)
@ -56,13 +59,24 @@ 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)
! Determine displacements for MPI_BANK type
call MPI_GET_ADDRESS(b % uid, bank_disp(1), ierr)
call MPI_GET_ADDRESS(b % xyz, bank_disp(2), ierr)
call MPI_GET_ADDRESS(b % uvw, bank_disp(3), ierr)
call MPI_GET_ADDRESS(b % E, bank_disp(4), ierr)
! Adjust displacements
base = bank_disp(1)
do i = 1, 4
bank_disp(i) = bank_disp(i) - base
end do
! Define MPI_BANK for fission sites
bank_blocks = (/ 1, 3, 3, 1 /)
bank_types = (/ MPI_INTEGER, 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)
end subroutine setup_mpi
@ -90,6 +104,10 @@ contains
& 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..."
call message(msg, 8)
! Determine starting index for fission bank and total sites in
! fission bank
@ -115,8 +133,11 @@ contains
! 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 there are less than n_particles particles banked,
! automatically add int(n_particles/total) sites to
! temp_sites. For example, if you need 1000 and 300 were
! banked, this would add 3 source sites per banked site and the
! remaining 100 would be randomly sampled.
if (total < n_particles) then
do j = 1,int(n_particles/total)
index = index + 1
@ -139,7 +160,7 @@ contains
temp_sites(count) = fission_bank(i-start)
temp_sites(count) % uid = index
end if
sites_needed = sites_needed-1
sites_needed = sites_needed - 1
end if
end do
@ -154,18 +175,18 @@ contains
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, &
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, &
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, &
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, &
call MPI_SEND(temp_sites(1), send_to_left, MPI_BANK, rank-1, 1, &
& MPI_COMM_WORLD, ierr)
end if
@ -192,6 +213,9 @@ contains
call copy_from_bank(temp_sites, i+1, j)
call copy_from_bank(right_bank, i+j+1, k)
end if
! Reset source index
source_index = 0
deallocate(left_bank )
deallocate(right_bank)

176
src/physics.f90
View file

@ -46,33 +46,34 @@ contains
end if
end if
if (verbosity >= 10) then
if (verbosity >= 9) then
msg = "=== Particle " // trim(int_to_str(p % uid)) // " ==="
call message(msg, 10)
end if
if (verbosity >= 10) then
i = cells(p % cell)%uid
msg = " Born in cell " // trim(int_to_str(i))
call message(msg, 10)
end if
! sample energy from Watt fission energy spectrum for U-235
p % E = watt_spectrum(0.988_8, 2.249_8)
! find energy index, interpolation factor
call find_energy_index(p)
IE = p % IE
f = p % interp
! Determine material total cross-section
Sigma = f*cMaterial%total_xs(IE) + (1-f)*cMaterial%total_xs(IE+1)
do while (p % alive)
! Copy energy index and interpolation factor
IE = p % IE
f = p % interp
! Determine material total cross-section
Sigma = (1-f)*cMaterial%total_xs(IE) + f*cMaterial%total_xs(IE+1)
! Find the distance to the nearest boundary
call dist_to_boundary(p, d_to_boundary, surf)
! Sample a distance to collision
d_to_collision = -log(rang()) / 1.0 ! Sigma
d_to_collision = -log(rang()) / Sigma
! Select smaller of the two distances
distance = min(d_to_boundary, d_to_collision)
@ -169,7 +170,7 @@ contains
IE = table%grid_index(p % IE)
f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE))
Sigma = density_i*(f*table%sigma_t(IE) + (1-f)*(table%sigma_t(IE+1)))
Sigma = density_i*((ONE-f)*table%Sigma_t(IE) + f*(table%Sigma_t(IE+1)))
Sigma_t(i) = Sigma
end do
@ -199,8 +200,8 @@ contains
rxn => table%reactions(i)
if (rxn%MT >= 200) cycle
if (IE < rxn%IE) cycle
prob = prob + density * (f*rxn%sigma(IE-rxn%IE+1) &
& + (1-f)*(rxn%sigma(IE-rxn%IE+2)))
prob = prob + density * ((ONE-f)*rxn%sigma(IE-rxn%IE+1) &
& + f*(rxn%sigma(IE-rxn%IE+2)))
if (r1 < prob) exit
end do
if (verbosity >= 10) then
@ -215,18 +216,19 @@ contains
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)
case (22, 23, 28, 29, 32:36, 43, 44, 51:91)
call inelastic_scatter(p, table, rxn)
case (18:21, 38)
call n_fission(p, table, rxn)
case (51:90)
call level_scatter(p, table, rxn)
case (102)
call n_gamma(p)
case default
call elastic_scatter(p, table%awr)
end select
! find energy index, interpolation factor
call find_energy_index(p)
end subroutine collision
!=====================================================================
@ -290,9 +292,6 @@ contains
p % uvw(2) = vy/vel
p % uvw(3) = vz/vel
! find energy index, interpolation factor
call find_energy_index(p)
end subroutine elastic_scatter
!=====================================================================
@ -365,10 +364,10 @@ contains
! appropriate index and interpolation factor
if (E < table % nu_t_data(loc+1)) then
j = 1
f = 0.0
f = ZERO
elseif (E > table % nu_t_data(loc+NE)) then
j = NE - 1
f = 1.0
f = ONE
else
j = binary_search(table % nu_t_data(loc+1), NE, E)
f = (E - table % nu_t_data(loc+j)) / &
@ -414,10 +413,10 @@ contains
! appropriate index and interpolation factor
if (E < table % nu_p_data(loc+1)) then
j = 1
f = 0.0
f = ZERO
elseif (E > table % nu_p_data(loc+NE)) then
j = NE - 1
f = 1.0
f = ONE
else
j = binary_search(table % nu_p_data(loc+1), NE, E)
f = (E - table % nu_p_data(loc+j)) / &
@ -451,10 +450,10 @@ contains
! appropriate index and interpolation factor
if (E < table % nu_d_data(loc+1)) then
j = 1
f = 0.0
f = ZERO
elseif (E > table % nu_d_data(loc+NE)) then
j = NE - 1
f = 1.0
f = ONE
else
j = binary_search(table % nu_d_data(loc+1), NE, E)
f = (E - table % nu_d_data(loc+j)) / &
@ -472,7 +471,7 @@ contains
! TODO: Heat generation from fission
! Sample number of neutrons produced
nu_total = p % wgt * nu_total
nu_total = p%wgt / keff * nu_total
if (rang() > nu_total - int(nu_total)) then
nu = int(nu_total)
else
@ -487,17 +486,13 @@ contains
fission_bank(i) % xyz = p % xyz
! sample cosine of angle
if (rxn % has_angle_dist) then
mu = sample_angle(rxn, E)
else
mu = TWO * rang() - ONE
end if
mu = sample_angle(rxn, E)
! Sample azimuthal angle uniformly in [0,2*pi)
phi = TWO*PI*rang()
fission_bank(i) % uvw(1) = mu
fission_bank(i) % uvw(2) = sqrt(1. - mu**2) * cos(phi)
fission_bank(i) % uvw(3) = sqrt(1. - mu**2) * sin(phi)
fission_bank(i) % uvw(2) = sqrt(ONE - mu**2) * cos(phi)
fission_bank(i) % uvw(3) = sqrt(ONE - mu**2) * sin(phi)
! determine energy of fission neutron
call sample_energy(rxn, E, E_out, mu_out)
@ -512,62 +507,6 @@ contains
end subroutine n_fission
!=====================================================================
! LEVEL_SCATTER
!=====================================================================
subroutine level_scatter(p, table, rxn)
type(Particle), pointer :: p
type(AceContinuous), pointer :: table
type(AceReaction), pointer :: rxn
real(8) :: A ! atomic weight ratio of nuclide
real(8) :: E_in ! incoming energy
real(8) :: mu_cm ! cosine of scattering angle in center-of-mass
real(8) :: mu_lab ! cosine of scattering angle in laboratory
real(8) :: E_cm ! outgoing energy in center-of-mass
real(8) :: E_lab ! outgoing energy in laboratory
character(250) :: msg ! error message
! copy energy of neutron
E_in = p % E
! determine A
A = table % awr
! determine if scattering is in CM (it should be!)
if (rxn % TY < 0) then
! scattering angle in center-of-mass
mu_cm = sample_angle(rxn, E_in)
else
msg = "Level inelastic scattering should not sample angle &
&in laboratory system!"
call error(msg)
end if
! sample outgoing energy in center-of-mass
call sample_energy(rxn, E_in, E_cm)
! determine outgoing energy in lab
E_lab = E_cm + (E_in + TWO * mu_cm * (A+ONE) * sqrt(E_in * E_cm)) &
& / ((A+ONE)*(A+ONE))
! determine outgoing angle in lab
mu_lab = mu_cm * sqrt(E_cm/E_lab) + ONE/(A+ONE) * sqrt(E_in/E_lab)
! change direction of particle
call rotate_angle(p, mu_lab)
! change energy of particle
p % E = E_lab
! find energy index, interpolation factor
call find_energy_index(p)
end subroutine level_scatter
!=====================================================================
! INELASTIC_SCATTER
!=====================================================================
@ -596,9 +535,7 @@ contains
law = rxn % edist % law
! sample scattering angle
if (rxn % has_angle_dist) then
mu = sample_angle(rxn, E_in)
end if
mu = sample_angle(rxn, E_in)
! sample outgoing energy
if (law == 44 .or. law == 61) then
@ -626,9 +563,6 @@ contains
! change energy of particle
p % E = E
! find energy index, interpolation factor
call find_energy_index(p)
end subroutine inelastic_scatter
!=====================================================================
@ -661,9 +595,7 @@ contains
law = rxn % edist % law
! sample scattering angle
if (rxn % has_angle_dist) then
mu = sample_angle(rxn, E_in)
end if
mu = sample_angle(rxn, E_in)
! sample outgoing energy
if (law == 44 .or. law == 61) then
@ -691,9 +623,6 @@ contains
! 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
@ -751,6 +680,13 @@ contains
real(8) :: p0,p1 ! probability distribution
character(250) :: msg ! error message
! check if reaction has angular distribution -- if not, sample
! outgoing angle isotropically
if (.not. rxn % has_angle_dist) then
mu = TWO * rang() - ONE
return
end if
! determine number of incoming energies
n = rxn % adist % n_energy
@ -759,10 +695,10 @@ contains
! the first or last bins
if (E < rxn % adist % energy(1)) then
i = 1
r = 0.0
r = ZERO
elseif (E > rxn % adist % energy(n)) then
i = n - 1
r = 1.0
r = ONE
else
i = binary_search(rxn % adist % energy, n, E)
r = (E - rxn % adist % energy(i)) / &
@ -1010,10 +946,10 @@ contains
loc = 2 + 2*NR
if (E_in < rxn % edist % data(loc+1)) then
i = 1
r = 0.0
r = ZERO
elseif (E_in > rxn % edist % data(loc+NE)) then
i = NE - 1
r = 1.0
r = ONE
else
i = binary_search(rxn % edist % data(loc+1), NE, E_in)
r = (E_in - rxn%edist%data(loc+i)) / &
@ -1123,10 +1059,10 @@ contains
loc = 2 + 2*NR
if (E_in < rxn % edist % data(loc+1)) then
i = 1
r = 0.0
r = ZERO
elseif (E_in > rxn % edist % data(loc+NE)) then
i = NE - 1
r = 1.0
r = ONE
else
i = binary_search(rxn % edist % data(loc+1), NE, E_in)
r = (E_in - rxn%edist%data(loc+i)) / &
@ -1160,10 +1096,10 @@ contains
loc = 2 + 2*NR
if (E_in < rxn % edist % data(loc+1)) then
i = 1
r = 0.0
r = ZERO
elseif (E_in > rxn % edist % data(loc+NE)) then
i = NE - 1
r = 1.0
r = ONE
else
i = binary_search(rxn % edist % data(loc+1), NE, E_in)
r = (E_in - rxn%edist%data(loc+i)) / &
@ -1175,6 +1111,10 @@ contains
T = rxn%edist%data(loc+i) + r * &
& (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i))
! determine restriction energy
loc = loc + NE
U = rxn % edist % data(loc + 1)
! sample outgoing energy based on evaporation spectrum
! probability density function
do
@ -1202,10 +1142,10 @@ contains
loc = 2 + 2*NR
if (E_in < rxn % edist % data(loc+1)) then
i = 1
r = 0.0
r = ZERO
elseif (E_in > rxn % edist % data(loc+NE)) then
i = NE - 1
r = 1.0
r = ONE
else
i = binary_search(rxn % edist % data(loc+1), NE, E_in)
r = (E_in - rxn%edist%data(loc+i)) / &
@ -1232,10 +1172,10 @@ contains
loc = loc + 2 + 2*NR
if (E_in < rxn % edist % data(loc+1)) then
i = 1
r = 0.0
r = ZERO
elseif (E_in > rxn % edist % data(loc+NE)) then
i = NE - 1
r = 1.0
r = ONE
else
i = binary_search(rxn % edist % data(loc+1), NE, E_in)
r = (E_in - rxn%edist%data(loc+i)) / &
@ -1274,10 +1214,10 @@ contains
loc = 2 + 2*NR
if (E_in < rxn % edist % data(loc+1)) then
i = 1
r = 0.0
r = ZERO
elseif (E_in > rxn % edist % data(loc+NE)) then
i = NE - 1
r = 1.0
r = ONE
else
i = binary_search(rxn % edist % data(loc+1), NE, E_in)
r = (E_in - rxn%edist%data(loc+i)) / &
@ -1419,10 +1359,10 @@ contains
loc = 2 + 2*NR
if (E_in < rxn % edist % data(loc+1)) then
i = 1
r = 0.0
r = ZERO
elseif (E_in > rxn % edist % data(loc+NE)) then
i = NE - 1
r = 1.0
r = ONE
else
i = binary_search(rxn % edist % data(loc+1), NE, E_in)
r = (E_in - rxn%edist%data(loc+i)) / &

4
src/source.f90
View file

@ -3,6 +3,7 @@ module source
use global
use mcnp_random, only: rang, RN_init_particle
use output, only: message
use physics, only: watt_spectrum
implicit none
@ -78,6 +79,9 @@ contains
p % universe = 0
p % wgt = ONE
p % alive = .true.
! sample energy from Watt fission energy spectrum for U-235
p % E = watt_spectrum(0.988_8, 2.249_8)
end do
end if
end do

68
src/tallies.f90 Normal file
View file

@ -0,0 +1,68 @@
module tallies
use global
use mpi
use output, only: message, error
implicit none
contains
!=====================================================================
! CALCULATE_KEFF
!=====================================================================
subroutine calculate_keff(i_cycle)
integer, intent(in) :: i_cycle ! index of current cycle
integer :: total_bank
integer :: n
integer :: ierr
real(8) :: kcoll ! keff collision estimator
real(8) :: ktemp ! MPI-reduced keff and stdev
real(8), save :: k1 = 0. ! accumulated keff
real(8), save :: k2 = 0. ! accumulated keff**2
real(8) :: std ! stdev of keff over active cycles
character(250) :: msg
msg = "Calculate cycle keff..."
call message(msg, 8)
! set k1 and k2 at beginning of run
if (i_cycle == 1) then
k1 = ZERO
k2 = ZERO
end if
! Collect number bank sites onto master process
call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER, MPI_SUM, 0, &
& MPI_COMM_WORLD, ierr)
! Collect statistics and print output
if (master) then
kcoll = real(total_bank)/real(n_particles)*keff
if (i_cycle > n_inactive) then
n = i_cycle - n_inactive
k1 = k1 + kcoll
k2 = k2 + kcoll**2
keff = k1/n
std = sqrt((k2/n-keff**2)/n)
if (i_cycle > n_inactive+1) then
write(6,101) i_cycle, kcoll, keff, std
else
write(6,100) i_cycle, kcoll
end if
else
write(6,100) i_cycle, kcoll
keff = kcoll
end if
end if
call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, ierr)
100 format (2X,I4,3X,F7.5)
101 format (2X,I4,3X,F7.5,10X,F7.5,2X,F7.5)
end subroutine calculate_keff
end module tallies

55
src/unittest.f90
View file

@ -1,9 +1,12 @@
program unittest
use global
use types, only: Bank
use energy_grid, only: add_grid_points
use mpi_routines, only: setup_mpi, MPI_BANK
use mpi
call test_energy_grid()
call test_mpi_bank()
contains
@ -14,8 +17,8 @@ contains
subroutine test_energy_grid()
type(LinkedListGrid), pointer :: list => null()
type(LinkedListGrid), pointer :: current => null()
type(ListReal), pointer :: list => null()
type(ListReal), pointer :: current => null()
real(8), allocatable :: energy(:)
integer :: n, i
integer :: n_list
@ -29,16 +32,16 @@ contains
! create list
allocate(list)
current => list
current%energy = 4.23_8
current%data = 4.23_8
allocate(current%next)
current => current%next
current%energy = 7.23_8
current%data = 7.23_8
call add_grid_points(list, energy, n)
current => list
do while(associated(current))
print *,current%energy
print *,current%data
current => current%next
end do
@ -46,4 +49,44 @@ contains
end subroutine test_energy_grid
!=====================================================================
! TEST_MPI_BANK is a test of the derived datatype MPI_BANK that is
! used to send fission bank sites during synchronization
!=====================================================================
subroutine test_mpi_bank()
integer :: ierr
type(Bank) :: temp_bank(2)
call setup_mpi()
if (rank == 0) then
temp_bank(1) % uid = 123
temp_bank(1) % xyz = (/ 13.0_8, 9.0_8, 8.0_8 /)
temp_bank(1) % uvw = (/ 7.5_8, 6.0_8, 2.0_8 /)
temp_bank(1) % E = 4.592834_8
temp_bank(2) % uid = 321
temp_bank(2) % xyz = (/ 23.0_8, 29.0_8, 28.0_8 /)
temp_bank(2) % uvw = (/ 1.0_8, 9.0_8, 11.0_8 /)
temp_bank(2) % E = 0.327438_8
call MPI_SEND(temp_bank, 2, MPI_BANK, 1, 3, MPI_COMM_WORLD, ierr)
elseif (rank == 1) then
call MPI_RECV(temp_bank, 2, MPI_BANK, 0, 3, MPI_COMM_WORLD, &
& MPI_STATUS_IGNORE, ierr)
print *, temp_bank(1) % uid
print *, temp_bank(1) % xyz
print *, temp_bank(1) % uvw
print *, temp_bank(1) % E
print *, temp_bank(2) % uid
print *, temp_bank(2) % xyz
print *, temp_bank(2) % uvw
print *, temp_bank(2) % E
end if
call MPI_FINALIZE(ierr)
end subroutine test_mpi_bank
end program unittest