From eee77ec2f3e82d99c597f5634b575ba9c6892b76 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 22 Mar 2011 14:55:26 +0000 Subject: [PATCH] Fixed major bugs, namely -- material total xs, MPI_BANK, evaporation spectrum, interpolation. Added initial keff ability. Added Jezebel benchmark. --- ChangeLog | 22 +++++ benchmarks/pu-met-fast-001 | 24 +++++ src/Makefile | 10 ++- src/cross_section.f90 | 57 ++++++++---- src/global.f90 | 3 + src/input_sample | 2 +- src/main.f90 | 10 ++- src/mpi_routines.f90 | 66 +++++++++----- src/physics.f90 | 176 ++++++++++++------------------------- src/source.f90 | 4 + src/tallies.f90 | 68 ++++++++++++++ src/unittest.f90 | 55 ++++++++++-- 12 files changed, 328 insertions(+), 169 deletions(-) create mode 100644 benchmarks/pu-met-fast-001 create mode 100644 src/tallies.f90 diff --git a/ChangeLog b/ChangeLog index 47de20a9cc..90b9029d36 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,25 @@ +2011-03-22 Paul Romano + + * 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 * ace.f90: Adjusted locators energy and angle distributions for diff --git a/benchmarks/pu-met-fast-001 b/benchmarks/pu-met-fast-001 new file mode 100644 index 0000000000..e14789b57c --- /dev/null +++ b/benchmarks/pu-met-fast-001 @@ -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 diff --git a/src/Makefile b/src/Makefile index cbf9558a18..b9e0b2e165 100644 --- a/src/Makefile +++ b/src/Makefile @@ -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 diff --git a/src/cross_section.f90 b/src/cross_section.f90 index 4eafe7cafc..dc2ce8c90d 100644 --- a/src/cross_section.f90 +++ b/src/cross_section.f90 @@ -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 diff --git a/src/global.f90 b/src/global.f90 index de856ed8e2..a1c3dec0b6 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -200,6 +200,9 @@ contains ! Default verbosity verbosity = 5 + ! Defualt multiplication factor + keff = ONE + end subroutine set_defaults !===================================================================== diff --git a/src/input_sample b/src/input_sample index df77b711f2..53c7b920f7 100644 --- a/src/input_sample +++ b/src/input_sample @@ -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 diff --git a/src/main.f90 b/src/main.f90 index 287d98ffaa..336433c236 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -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 diff --git a/src/mpi_routines.f90 b/src/mpi_routines.f90 index 350ba0e860..bcf3c77695 100644 --- a/src/mpi_routines.f90 +++ b/src/mpi_routines.f90 @@ -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) diff --git a/src/physics.f90 b/src/physics.f90 index 47aa44a1f0..b7efa242e6 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -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)) / & diff --git a/src/source.f90 b/src/source.f90 index e1d83208d0..783fb03d55 100644 --- a/src/source.f90 +++ b/src/source.f90 @@ -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 diff --git a/src/tallies.f90 b/src/tallies.f90 new file mode 100644 index 0000000000..4b73322e8e --- /dev/null +++ b/src/tallies.f90 @@ -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 diff --git a/src/unittest.f90 b/src/unittest.f90 index 07c4426c89..afa46cefa2 100644 --- a/src/unittest.f90 +++ b/src/unittest.f90 @@ -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