diff --git a/ChangeLog b/ChangeLog index 211062cecd..1842b7db79 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,24 @@ +2011-04-25 Paul Romano + + * ace.f90: Fixed indentation of get_int and get_real. Added + subroutine get_macro_xs. + * fileio.f90: Added ability to read tally cards through read_tally + subroutine. Changed ERROR_CODE to ERROR_INT. + * global.f90: Added variables for tallies. Changed ERROR_CODE to + ERROR_INT and added ERROR_REAL. + * input_sample: Added tally card. + * main.f90: Use tallies_on variable when appropriate. + * output.f90: Add preliminary tally printing capabilites through + print_tally subroutine. + * physics.f90: Add call to score_tally for tallies. Changed + Sigma_t to Sigam_rxn. + * score.f90: Added score_tally and add_to_score subroutine. + * search.f90: Added interpolate subroutine. + * string.f90: Added is_number function. + * tallies.f90: CHANGED NAME TO score.f90 + * types.f90: Added list of tallies to Universe, Lattice and + Cell. Modified Tally type. + 2011-04-22 Paul Romano * geometry.f90: Deallocated count lists in neighbor_lists. diff --git a/src/Makefile b/src/Makefile index c6d7c828ab..4593ebe290 100644 --- a/src/Makefile +++ b/src/Makefile @@ -17,7 +17,7 @@ modules = ace.f90 \ search.f90 \ source.f90 \ string.f90 \ - tallies.f90 \ + score.f90 \ types.f90 main_objects = $(modules:.f90=.o) $(main:.f90=.o) @@ -70,14 +70,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 tallies.o + ace.o energy_grid.o score.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 endf.o score.o +score.o: global.o output.o ace.o search.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 mpi_routines.o diff --git a/src/ace.f90 b/src/ace.f90 index 1d90cd2619..39a9226746 100644 --- a/src/ace.f90 +++ b/src/ace.f90 @@ -1075,29 +1075,96 @@ contains ! in the XSS array !===================================================================== - function get_int(n_values) result(array) + function get_int(n_values) result(array) - integer, intent(in) :: n_values ! number of values to read - integer :: array(n_values) ! array of values + integer, intent(in) :: n_values ! number of values to read + integer :: array(n_values) ! array of values - array = int(XSS(XSS_index:XSS_index + n_values - 1)) - XSS_index = XSS_index + n_values - - end function get_int + array = int(XSS(XSS_index:XSS_index + n_values - 1)) + XSS_index = XSS_index + n_values + + end function get_int !===================================================================== ! GET_REAL returns an array of real(8)s read from the current position ! in the XSS array !===================================================================== - function get_real(n_values) result(array) + function get_real(n_values) result(array) - integer, intent(in) :: n_values ! number of values to read - real(8) :: array(n_values) ! array of values + integer, intent(in) :: n_values ! number of values to read + real(8) :: array(n_values) ! array of values - array = XSS(XSS_index:XSS_index + n_values - 1) - XSS_index = XSS_index + n_values + array = XSS(XSS_index:XSS_index + n_values - 1) + XSS_index = XSS_index + n_values - end function get_real + end function get_real + +!===================================================================== +! GET_MACRO_XS +!===================================================================== + + function get_macro_xs(p, mat, MT) result(xs) + + type(Particle), pointer :: p + type(Material), pointer :: mat + integer, intent(in) :: MT + real(8) :: xs + + integer :: i + integer :: n_isotopes + integer :: IE + real(8) :: density + real(8) :: density_i + real(8) :: sigma_i + real(8) :: f + type(AceContinuous), pointer :: table => null() + type(AceReaction), pointer :: rxn => null() + + ! initialize xs + xs = ZERO + + ! find material atom density + density = cMaterial%atom_density + + ! loop over all isotopes in material + n_isotopes = cMaterial % n_isotopes + do i = 1, n_isotopes + table => xs_continuous(cMaterial % table(i)) + + ! determine nuclide atom density + density_i = cMaterial%atom_percent(i) * density + + ! search nuclide energy grid + IE = table%grid_index(p % IE) + f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE)) + + ! handle special case of total cross section + if (MT == 1) then + xs = xs + density * (ONE-f) * table%sigma_t(IE) + & + & f * (table%sigma_t(IE+1)) + cycle + end if + + ! loop over reactions in isotope + do j = 1, table % n_reaction + rxn => table % reactions(i) + + ! check for matching MT + if (MT /= rxn % MT) cycle + + ! if energy is below threshold for this reaction, skip it + if (IE < rxn % IE) cycle + + ! add to cumulative probability + sigma_i = (ONE-f) * rxn%sigma(IE-rxn%IE+1) + & + & f * (rxn%sigma(IE-rxn%IE+2)) + end do + + ! calculate nuclide macroscopic cross-section + xs = xs + density_i * sigma_i + end do + + end function get_macro_xs end module ace diff --git a/src/fileio.f90 b/src/fileio.f90 index 232ce1b9bb..cd33598620 100644 --- a/src/fileio.f90 +++ b/src/fileio.f90 @@ -2,7 +2,8 @@ module fileio use global use types, only: Cell, Surface, ExtSource, ListKeyValueCI - use string, only: split_string_wl, lower_case, split_string, concatenate + use string, only: split_string_wl, lower_case, split_string, & + & concatenate, is_number use output, only: message, warning, error use data_structures, only: dict_create, dict_add_key, dict_get_key, & & dict_has_key, DICT_NULL, dict_keys, list_size @@ -12,7 +13,6 @@ module fileio type(DictionaryII), pointer :: & ! used to count how many cells each & ucount_dict => null(), & ! universe contains & bc_dict => null() ! store boundary conditions - integer, allocatable :: index_cell_in_univ(:) @@ -93,6 +93,7 @@ contains n_universes = 0 n_surfaces = 0 n_materials = 0 + n_tallies = 0 call dict_create(ucount_dict) call dict_create(universe_dict) @@ -131,6 +132,8 @@ contains n_surfaces = n_surfaces + 1 case ('material') n_materials = n_materials + 1 + case ('tally') + n_tallies = n_tallies + 1 case ('lattice') n_lattices = n_lattices + 1 @@ -167,6 +170,7 @@ contains allocate(lattices(n_lattices)) allocate(surfaces(n_surfaces)) allocate(materials(n_materials)) + allocate(tallies(n_tallies)) ! Also allocate a list for keeping track of where cells have been ! assigned in each universe @@ -178,6 +182,7 @@ contains call dict_create(surface_dict) call dict_create(bc_dict) call dict_create(material_dict) + call dict_create(tally_dict) ! We also need to allocate the cell count lists for each ! universe. The logic for this is a little more convoluted. In @@ -227,6 +232,7 @@ contains integer :: index_lattice ! index in lattices array integer :: index_material ! index in materials array integer :: index_source ! index in source array (?) + integer :: index_tally ! index in tally array character(250) :: line ! a line of words character(250) :: msg ! output/error message character(32) :: words(max_words) ! words on a single line @@ -239,6 +245,7 @@ contains index_surface = 0 index_lattice = 0 index_material = 0 + index_tally = 0 index_source = 0 open(FILE=filename, UNIT=in, STATUS='old', & @@ -280,6 +287,11 @@ contains index_material = index_material + 1 call read_material(index_material, words, n) + case ('tally') + ! Read tally entry + index_tally = index_tally + 1 + call read_tally(index_tally, words, n) + case ('source') call read_source(words, n) @@ -291,7 +303,7 @@ contains case ('verbosity') verbosity = str_to_int(words(2)) - if (verbosity == ERROR_CODE) then + if (verbosity == ERROR_INT) then msg = "Invalid verbosity: " // words(2) call error(msg) end if @@ -468,7 +480,7 @@ contains ! Read cell identifier c % uid = str_to_int(words(2)) - if (c % uid == ERROR_CODE) then + if (c % uid == ERROR_INT) then msg = "Invalid cell name: " // words(2) call error(msg) end if @@ -476,7 +488,7 @@ contains ! Read cell universe universe_num = str_to_int(words(3)) - if (universe_num == ERROR_CODE) then + if (universe_num == ERROR_INT) then msg = "Invalid universe: " // words(3) call error(msg) end if @@ -489,7 +501,7 @@ contains ! find universe universe_num = str_to_int(words(5)) - if (universe_num == ERROR_CODE) then + if (universe_num == ERROR_INT) then msg = "Invalid universe fill: " // words(5) call error(msg) end if @@ -506,7 +518,7 @@ contains c % type = CELL_NORMAL c % material = str_to_int(words(4)) c % fill = 0 - if (c % material == ERROR_CODE) then + if (c % material == ERROR_INT) then msg = "Invalid material number: " // words(4) call error(msg) end if @@ -647,7 +659,7 @@ contains ! Read surface identifier surface_uid = str_to_int(words(2)) - if (surface_uid == ERROR_CODE) then + if (surface_uid == ERROR_INT) then msg = "Invalid surface name: " // words(2) call error(msg) end if @@ -695,7 +707,7 @@ contains ! Read lattice universe universe_num = str_to_int(words(2)) - if (universe_num == ERROR_CODE) then + if (universe_num == ERROR_INT) then msg = "Invalid universe: " // words(2) call error(msg) end if @@ -717,10 +729,10 @@ contains ! Read number of lattice cells in each direction n_x = str_to_int(words(4)) n_y = str_to_int(words(5)) - if (n_x == ERROR_CODE) then + if (n_x == ERROR_INT) then msg = "Invalid number of lattice cells in x-direction: " // words(4) call error(msg) - elseif (n_y == ERROR_CODE) then + elseif (n_y == ERROR_INT) then msg = "Invalid number of lattice cells in y-direction: " // words(5) call error(msg) end if @@ -748,7 +760,7 @@ contains do i = 1, n_x index_word = 9 + j*n_x + i universe_num = str_to_int(words(index_word)) - if (universe_num == ERROR_CODE) then + if (universe_num == ERROR_INT) then msg = "Invalid universe number: " // words(index_word) call error(msg) end if @@ -799,6 +811,161 @@ contains end subroutine read_source +!===================================================================== +! READ_TALLY +!===================================================================== + + subroutine read_tally(index, words, n_words) + + integer, intent(in) :: index ! index in materials array + character(*), intent(in) :: words(n_words) ! words on material entry + integer, intent(in) :: n_words ! number of words + + integer :: i ! index in words array + integer :: j + integer :: k + integer :: count + integer :: MT + integer :: cell_uid + integer :: r_bins, c_bins, e_bins + real(8) :: E + character(32) :: word + character(250) :: msg + type(Tally), pointer :: t => null() + + t => tallies(index) + + ! Read tally identifier + t % uid = str_to_int(words(2)) + if (t % uid == ERROR_INT) then + msg = "Invalid tally name: " // words(2) + call error(msg) + end if + call dict_add_key(tally_dict, t % uid, index) + + i = 3 + do while (i <= n_words) + word = words(i) + call lower_case(word) + + select case (trim(word)) + case ('reaction') + ! ========================================================== + ! READ REACTION LIST + + ! Determine how many reactions are listed + do j = i+1, n_words + if (.not. is_number(words(j))) exit + end do + count = j - (i+1) + + ! Allocate space to store reactions + allocate(t % reactions(count)) + + ! Read reaction MT values + do j = 1, count + MT = str_to_int(words(i+j)) + if (MT == ERROR_INT) then + msg = "Invalid reaction MT on tally: " // int_to_str(t%uid) + call error(msg) + end if + t % reactions(j) = MT + end do + + ! Determine whether to sum reactions or individually bin + if (trim(words(i+count+1)) == 'sum') then + t % reaction_type = TALLY_SUM + r_bins = 1 + else + t % reaction_type = TALLY_BINS + r_bins = count + end if + + case ('cell') + ! ========================================================== + ! READ CELL LIST + + ! Determine how many reactions are listed + do j = i+1, n_words + if (.not. is_number(words(j))) exit + end do + count = j - (i+1) + + ! Allocate space to store reactions + allocate(t % cells(count)) + + ! Read reaction MT values + do j = 1, count + cell_uid = str_to_int(words(i+j)) + if (cell_uid == ERROR_INT) then + msg = "Invalid cell number on tally: " // int_to_str(t%uid) + call error(msg) + end if + t % cells(j) = cell_uid + end do + + ! Determine whether to sum reactions or individually bin + if (trim(words(i+count+1)) == 'sum') then + t % cell_type = TALLY_SUM + c_bins = 1 + else + t % cell_type = TALLY_BINS + c_bins = count + end if + + case ('energy') + ! ========================================================== + ! READ ENERGY LIST + + ! Determine how many energies are listed + do j = i+1, n_words + if (.not. is_number(words(j)(1:1))) exit + end do + count = j - (i+1) + + ! Allocate space to store energies + allocate(t % energies(count)) + + ! Read reaction MT values + do j = 1, count + E = str_to_real(words(i+j)) + if (E == ERROR_REAL) then + msg = "Invalid energy on tally: " // int_to_str(t % uid) + call error(msg) + end if + t % energies(j) = E + end do + + e_bins = count - 1 + + case ('material') + ! TODO: Add ability to read material lists + + case ('universe') + ! TODO: Add ability to read universe lists + + case ('lattice') + ! TODO: Add lattice tally ability + + end select + + ! Move index in words forward + i = i + 1 + count + + end do + + ! allocate tally scores + allocate(t % score(r_bins, c_bins, e_bins)) + + ! initialize tallies + forall (i=1:r_bins, j=1:c_bins, k=1:e_bins) + t % score(i,j,k) % n_events = 0 + t % score(i,j,k) % val = ZERO + t % score(i,j,k) % val_sq = ZERO + end forall + + end subroutine read_tally + !===================================================================== ! READ_MATERIAL parses a material card. Note that atom percents and ! densities are normalized in a separate routine @@ -964,21 +1131,21 @@ contains ! Read number of cycles n_cycles = str_to_int(words(2)) - if (n_cycles == ERROR_CODE) then + if (n_cycles == ERROR_INT) then msg = "Invalid number of cycles: " // words(2) call error(msg) end if ! Read number of inactive cycles n_inactive = str_to_int(words(3)) - if (n_inactive == ERROR_CODE) then + if (n_inactive == ERROR_INT) then msg = "Invalid number of inactive cycles: " // words(2) call error(msg) end if ! Read number of particles n_particles = str_to_int(words(4)) - if (n_particles == ERROR_CODE) then + if (n_particles == ERROR_INT) then msg = "Invalid number of particles: " // words(2) call error(msg) end if diff --git a/src/global.f90 b/src/global.f90 index b19e72730e..48852565d2 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -16,11 +16,15 @@ module global type(Material), allocatable, target :: materials(:) type(Isotope), allocatable, target :: isotopes(:) type(xsData), allocatable, target :: xsdatas(:) - integer :: n_cells ! # of cells - integer :: n_universes ! # of universes - integer :: n_lattices ! # of lattices - integer :: n_surfaces ! # of surfaces - integer :: n_materials ! # of materials + type(Tally), allocatable, target :: tallies(:) + type(Tally), allocatable, target :: tallies_global(:) + integer :: n_cells ! # of cells + integer :: n_universes ! # of universes + integer :: n_lattices ! # of lattices + integer :: n_surfaces ! # of surfaces + integer :: n_materials ! # of materials + integer :: n_tallies ! # of tallies + integer :: n_tallies_global ! # of global tallies ! These dictionaries provide a fast lookup mechanism type(DictionaryII), pointer :: cell_dict @@ -29,6 +33,7 @@ module global type(DictionaryII), pointer :: surface_dict type(DictionaryII), pointer :: material_dict type(DictionaryII), pointer :: isotope_dict + type(DictionaryII), pointer :: tally_dict type(DictionaryCI), pointer :: xsdata_dict type(DictionaryCI), pointer :: ace_dict @@ -68,6 +73,8 @@ module global ! cycle keff real(8) :: keff + logical :: tallies_on + ! Parallel processing variables integer :: n_procs ! number of processes integer :: rank ! rank of process @@ -223,6 +230,13 @@ module global & N_PT = 116, & & N_DA = 117 + ! Tally distinguishers + integer, parameter :: & + & TALLY_FLUX = -1, & ! tally the flux only + & TALLY_ALL = 0, & ! tally all reactions + & TALLY_BINS = 1, & ! tally individual (reactions, cells) + & TALLY_SUM = 2 ! tally sum of specified (reactions, cells) + ! Fission neutron emission (nu) type integer, parameter :: & & NU_NONE = 0, & ! No nu values (non-fissionable) @@ -232,7 +246,10 @@ module global ! Integer code for read error -- better hope this number is never ! used in an input file! integer, parameter :: & - & ERROR_CODE = -huge(0) + & ERROR_INT = -huge(0) + + real(8), parameter :: & + & ERROR_REAL = -huge(0.0_8) * 0.917826354_8 ! The verbosity controls how much information will be printed to the ! screen and in logs @@ -358,7 +375,7 @@ contains ! read string into integer read(UNIT=str, FMT=fmt, IOSTAT=ioError) num - if (ioError > 0) num = ERROR_CODE + if (ioError > 0) num = ERROR_INT end function str_to_int @@ -409,9 +426,7 @@ contains ! Read string read(string, fmt, iostat=readError) num - if (readError > 0) then - ! return error code for real? - end if + if (readError > 0) num = ERROR_REAL end function str_to_real diff --git a/src/input_sample b/src/input_sample index 8d037878dd..575e0d7338 100644 --- a/src/input_sample +++ b/src/input_sample @@ -20,6 +20,8 @@ material 41 -1.0 & source box -4 -4 -4 4 4 4 +tally 1 reaction 2 cell 2 energy limits + xs_data /opt/serpent/xsdata/endfb7 criticality 15 5 10000 verbosity 7 diff --git a/src/main.f90 b/src/main.f90 index 5b6e35f743..51f9be5c87 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -13,7 +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, t_sync - use tallies, only: calculate_keff + use score, only: calculate_keff #ifdef MPI use mpi @@ -111,6 +111,8 @@ contains msg = "Running problem..." call message(msg, 6) + tallies_on = .false. + #ifdef MPI t0 = MPI_WTIME() #endif @@ -148,9 +150,10 @@ contains ! print cycle information - end do CYCLE_LOOP + ! Turn tallies on once inactive cycles are complete + if (i_cycle == n_inactive) tallies_on = .true. - ! Collect all tallies and print + end do CYCLE_LOOP #ifdef MPI ! print run time diff --git a/src/output.f90 b/src/output.f90 index 3fb3507979..4066f3daf4 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -511,6 +511,66 @@ contains end subroutine print_material +!===================================================================== +! PRINT_TALLY displays the attributes of a tally +!===================================================================== + + subroutine print_tally(tal) + + type(Tally), pointer :: tal + + integer :: i + integer :: MT + character(250) :: string + + write(ou,*) 'Tally ' // int_to_str(tal % uid) + + select case (tal % reaction_type) + case (TALLY_FLUX) + write(ou,*) ' Type: Flux' + case (TALLY_ALL) + write(ou,*) ' Type: Total Collision Rate' + case (TALLY_BINS) + write(ou,*) ' Type: Partial reactions' + case (TALLY_SUM) + write(ou,*) ' Type: Partial reactions (summed)' + end select + + select case (tal % cell_type) + case (TALLY_BINS) + write(ou,*) ' Cell Type: Separate bins' + case (TALLY_SUM) + write(ou,*) ' Cell Type: Sum over cells' + end select + + if (allocated(tal % reactions)) then + string = "" + do i = 1, size(tal % reactions) + MT = tal % reactions(i) + string = trim(string) // ' ' // trim(reaction_name(MT)) + end do + write(ou,*) ' Reactions:' // trim(string) + end if + + if (allocated(cells)) then + string = "" + do i = 1, size(tal % cells) + string = trim(string) // ' ' // trim(int_to_str(tal % cells(i))) + end do + write(ou,*) ' Cells:' // trim(string) + end if + + if (allocated(tal % energies)) then + string = "" + do i = 1, size(tal % energies) + string = trim(string) // ' ' // trim(real_to_str(tal % energies(i))) + end do + write(ou,*) ' Energies:' // trim(string) + end if + write(ou,*) + + end subroutine print_tally + !===================================================================== ! PRINT_SUMMARY displays the attributes of all cells, universes, ! surfaces and materials read in the input file. Very useful for @@ -524,6 +584,7 @@ contains type(Universe), pointer :: u => null() type(Lattice), pointer :: l => null() type(Material), pointer :: m => null() + type(Tally), pointer :: t => null() integer :: i ! print summary of cells @@ -578,11 +639,24 @@ contains call print_material(m) end do + ! print summary of tallies + if (n_tallies > 0) then + write(ou,*) '=============================================' + write(ou,*) '=> TALLY SUMMARY <=' + write(ou,*) '=============================================' + write(ou,*) + do i = 1, n_tallies + t=> tallies(i) + call print_tally(t) + end do + end if + nullify(s) nullify(c) nullify(u) nullify(l) nullify(m) + nullify(t) end subroutine print_summary diff --git a/src/physics.f90 b/src/physics.f90 index 25fcbcfd04..a2a0c0994a 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -8,6 +8,7 @@ module physics use output, only: error, warning, message, print_particle use search, only: binary_search use endf, only: reaction_name + use score, only: score_tally implicit none @@ -81,8 +82,6 @@ contains p%xyz = p%xyz + distance * p%uvw p%xyz_local = p%xyz_local + distance * p%uvw - ! Add pathlength tallies - if (d_to_collision > d_to_boundary) then p % cell = 0 if (in_lattice) then @@ -159,50 +158,51 @@ contains real(8) :: density_i ! atom density of nuclide real(8) :: prob ! cumulative probability real(8) :: r1 ! random number - real(8), allocatable :: Sigma_t(:) ! macroscopic xs for each nuclide + real(8) :: flux ! collision estimator of flux + real(8), allocatable :: Sigma_rxn(:) ! macroscopic xs for each nuclide character(250) :: msg ! output/error message type(AceContinuous), pointer :: table type(AceReaction), pointer :: rxn - ! TODO: tallies - density = cMaterial%atom_density ! calculate total cross-section for each nuclide at current energy ! in order to create discrete pdf for sampling nuclide n_isotopes = cMaterial%n_isotopes - allocate(Sigma_t(n_isotopes)) + allocate(Sigma_rxn(n_isotopes)) do i = 1, n_isotopes table => xs_continuous(cMaterial%table(i)) - density_i = cMaterial%atom_percent(i)*density + + ! determine nuclide atom density + density_i = cMaterial%atom_percent(i) * density ! search nuclide energy grid IE = table%grid_index(p % IE) f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE)) - ! calculate macroscopic xs for this nuclide - Sigma_t(i) = density_i*((ONE-f)*table%Sigma_t(IE) + & + ! calculate nuclide macroscopic cross-section + Sigma_rxn(i) = density_i*((ONE-f)*table%Sigma_t(IE) + & & f*(table%Sigma_t(IE+1))) end do - total = sum(Sigma_t) + total = sum(Sigma_rxn) ! sample nuclide r1 = rang() prob = ZERO do i = 1, n_isotopes - prob = prob + Sigma_t(i) / total + prob = prob + Sigma_rxn(i) / total if (r1 < prob) exit end do ! Get table, total xs, interpolation factor density_i = cMaterial%atom_percent(i)*density table => xs_continuous(cMaterial%table(i)) - sigma = Sigma_t(i) / density_i + sigma = Sigma_rxn(i) / density_i IE = table%grid_index(p % IE) f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE)) ! free memory - deallocate(Sigma_t) + deallocate(Sigma_rxn) ! sample reaction channel r1 = rang()*sigma @@ -248,6 +248,12 @@ contains call warning(msg) end select + ! Add collision estimator tallies + if (tallies_on) then + flux = p%wgt / SIGMA + call score_tally(p, flux) + end if + ! check for very low energy if (p % E < 1.0e-100_8) then p % alive = .false. diff --git a/src/score.f90 b/src/score.f90 new file mode 100644 index 0000000000..c1e9c67b7a --- /dev/null +++ b/src/score.f90 @@ -0,0 +1,277 @@ +module score + + use global + use output, only: message, error + use ace, only: get_macro_xs + use search, only: binary_search + +#ifdef MPI + use mpi +#endif + + implicit none + +contains + +!===================================================================== +! CALCULATE_KEFF +!===================================================================== + + subroutine calculate_keff(i_cycle) + + integer, intent(in) :: i_cycle ! index of current cycle + + integer(8) :: 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 + +#ifdef MPI + ! Collect number bank sites onto master process + call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & + & MPI_COMM_WORLD, ierr) +#else + total_bank = n_bank +#endif + + ! 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 + +#ifdef MPI + call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, ierr) +#endif + +100 format (2X,I4,3X,F7.5) +101 format (2X,I4,3X,F7.5,10X,F7.5,2X,F7.5) + + end subroutine calculate_keff + +!===================================================================== +! SCORE_TALLY +!===================================================================== + + subroutine score_tally(p, flux) + + type(Particle), pointer :: p ! particle + real(8), intent(in) :: flux ! estimator of flux + + integer :: i ! index for tallies in cell + integer :: j ! index for reactions in tally + integer :: e_bin ! energy bin + integer :: c_bin ! cell bin + integer :: r_bin ! reaction bin + integer :: n_energy ! number of energy bins + integer :: n_reaction ! number of reactions + integer :: type ! cell or reaction type + integer :: MT ! reaction MT value + real(8) :: E ! energy of particle + real(8) :: val ! value to score + real(8) :: Sigma ! macroscopic cross section of reaction + character(250) :: msg ! output/error message + type(Cell), pointer :: c => null() + type(Tally), pointer :: t => null() + + ! ================================================================ + ! HANDLE LOCAL TALLIES + + c => cells(p % cell) + if (.not. allocated(c % tallies)) return + + ! if so, loop over each tally + do i = 1, size(c % tallies) + t => tallies(c % tallies(i)) + + ! ============================================================= + ! DETERMINE ENERGY BIN + if (allocated(t % energies)) then + E = p % E + n_energy = size(t % energies) + if (E < t % energies(1) .or. E > t % energies(n_energy)) then + ! Energy outside of specified range + cycle + else + e_bin = binary_search(t % energies, n_energy, E) + end if + end if + + ! ============================================================= + ! DETERMINE CELL BIN + type = t % cell_type + if (type == TALLY_SUM) then + ! Sum tallies from separate cells into one bin + c_bin = 1 + elseif (type == TALLY_BINS) then + ! Need to determine cell bin + do j = 1, size(t % cells) + if (p % cell == t % cells(j)) then + c_bin = j + exit + end if + end do + else + msg = "Invalid type for cell bins in tally " // int_to_str(t % uid) + call error(msg) + end if + + ! ============================================================= + ! DETERMINE REACTION BIN AND ADD VALUE TO SCORE + type = t % reaction_type + if (type == TALLY_FLUX) then + ! Tally flux only + val = flux + call add_to_score(t % score(r_bin, c_bin, e_bin), val) + + elseif (type == TALLY_ALL) then + ! Tally total reaction rate + val = ONE + r_bin = 1 + call add_to_score(t % score(r_bin, c_bin, e_bin), val) + + elseif (type == TALLY_BINS) then + ! Individually bin reactions + n_reaction = t % reaction_type + + r_bin = 0 + do j = 1, n_reaction + MT = t % reactions(j) + Sigma = get_macro_xs(p, cMaterial, MT) + val = Sigma * flux + r_bin = r_bin + 1 + call add_to_score(t % score(r_bin, c_bin, e_bin), & + & val) + end do + elseif (type == TALLY_SUM) then + ! Tally reactions in one bin + n_reaction = t % reaction_type + + r_bin = 1 + do j = 1, n_reaction + MT = t % reactions(j) + Sigma = get_macro_xs(p, cMaterial, MT) + val = Sigma * flux + call add_to_score(t % score(r_bin, c_bin, e_bin), & + & val) + end do + end if + end do + + ! ================================================================ + ! HANDLE GLOBAL TALLIES + + do i = 1, n_tallies_global + t => tallies_global(i) + + ! ============================================================= + ! DETERMINE ENERGY BIN + if (allocated(t % energies)) then + E = p % E + n_energy = size(t % energies) + if (E < t % energies(1) .or. E > t % energies(n_energy)) then + ! Energy outside of specified range + cycle + else + e_bin = binary_search(t % energies, n_energy, E) + end if + end if + + ! Since it's a global tally, the cell bin is unity + c_bin = 1 + + ! ============================================================= + ! DETERMINE REACTION BIN AND ADD VALUE TO SCORE + type = t % reaction_type + if (type == TALLY_FLUX) then + ! Tally flux only + val = flux + call add_to_score(t % score(r_bin, c_bin, e_bin), val) + + elseif (type == TALLY_ALL) then + ! Tally total reaction rate + val = ONE + r_bin = 1 + call add_to_score(t % score(r_bin, c_bin, e_bin), val) + + elseif (type == TALLY_BINS) then + ! Individually bin reactions + n_reaction = t % reaction_type + + r_bin = 0 + do j = 1, n_reaction + MT = t % reactions(j) + Sigma = get_macro_xs(p, cMaterial, MT) + val = Sigma * flux + r_bin = r_bin + 1 + call add_to_score(t % score(r_bin, c_bin, e_bin), & + & val) + end do + elseif (type == TALLY_SUM) then + ! Tally reactions in one bin + n_reaction = t % reaction_type + + r_bin = 1 + do j = 1, n_reaction + MT = t % reactions(j) + Sigma = get_macro_xs(p, cMaterial, MT) + val = Sigma * flux + call add_to_score(t % score(r_bin, c_bin, e_bin), & + & val) + end do + end if + + end do + + ! ================================================================ + ! TODO: Add lattice tallies + + ! ================================================================ + ! TODO: Add mesh tallies + + end subroutine score_tally + +!===================================================================== +! ADD_TO_SCORE +!===================================================================== + + subroutine add_to_score(score, val) + + type(TallyScore), intent(inout) :: score + real(8), intent(in) :: val + + score % n_events = score % n_events + 1 + score % val = score % val + val + score % val_sq = score % val_sq + val*val + + end subroutine add_to_score + +end module score diff --git a/src/search.f90 b/src/search.f90 index 48e9b4b8fa..103808a54d 100644 --- a/src/search.f90 +++ b/src/search.f90 @@ -55,4 +55,20 @@ contains end function binary_search +!===================================================================== +! INTERPOLATE +!===================================================================== + + function interpolate(array, n, index, f) result(val) + + real(8), intent(in) :: array(n) + integer, intent(in) :: n + integer, intent(in) :: index + real(8), intent(in) :: f + real(8) :: val + + val = (ONE-f) * array(index) + f * array(index+1) + + end function interpolate + end module search diff --git a/src/string.f90 b/src/string.f90 index 94479a2059..03d67cb43f 100644 --- a/src/string.f90 +++ b/src/string.f90 @@ -150,7 +150,6 @@ contains !===================================================================== elemental subroutine lower_case(word) - ! convert a word to lower case character(*), intent(inout) :: word @@ -164,4 +163,25 @@ contains end subroutine lower_case +!===================================================================== +! IS_NUMBER determines whether a string of characters is all 0-9 +! characters +!===================================================================== + + function is_number(word) result(number) + + character(*), intent(in) :: word + logical :: number + + integer :: i + integer :: ic + + number = .true. + do i = 1, len_trim(word) + ic = ichar(word(i:i)) + if (ic < 48 .or. ic >= 58) number = .false. + end do + + end function is_number + end module string diff --git a/src/tallies.f90 b/src/tallies.f90 deleted file mode 100644 index be3d6ff497..0000000000 --- a/src/tallies.f90 +++ /dev/null @@ -1,78 +0,0 @@ -module tallies - - use global - use output, only: message, error - -#ifdef MPI - use mpi -#endif - - implicit none - -contains - -!===================================================================== -! CALCULATE_KEFF -!===================================================================== - - subroutine calculate_keff(i_cycle) - - integer, intent(in) :: i_cycle ! index of current cycle - - integer(8) :: 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 - -#ifdef MPI - ! Collect number bank sites onto master process - call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & - & MPI_COMM_WORLD, ierr) -#else - total_bank = n_bank -#endif - - ! 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 - -#ifdef MPI - call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, ierr) -#endif - -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/types.f90 b/src/types.f90 index ebfcb4c730..aef5c3b39d 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -15,6 +15,7 @@ module types real(8) :: x0 ! Translation in x-coordinate real(8) :: y0 ! Translation in y-coordinate real(8) :: z0 ! Translation in z-coordinate + integer, allocatable :: tallies(:) end type Universe !===================================================================== @@ -33,6 +34,7 @@ module types real(8) :: width_x ! width of lattice cell real(8) :: width_y ! width of lattice cell integer, allocatable :: element(:,:) ! specified universes + integer, allocatable :: tallies(:) end type Lattice !===================================================================== @@ -66,6 +68,7 @@ module types & surfaces(:) ! List of surfaces bounding cell -- note ! that parentheses, union, etc operators ! will be listed here too + integer, allocatable :: tallies(:) end type Cell !===================================================================== @@ -124,6 +127,9 @@ module types integer :: uid integer :: type real(8) :: volume + integer :: cell_type + integer :: reaction_type + integer :: material_type integer, allocatable :: reactions(:) integer, allocatable :: cells(:) integer, allocatable :: materials(:) @@ -134,7 +140,7 @@ module types integer :: n_x integer :: n_y integer :: n_z - type(TallyScore), allocatable :: score(:) + type(TallyScore), allocatable :: score(:,:,:) end type Tally !=====================================================================