From 03a2634fb4e1e2865745ec84d1e38a7dbac252c5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 2 Mar 2011 15:23:55 +0000 Subject: [PATCH] Added initial universe capability (not complete), and fixed angular distribution sampling. --- ChangeLog | 23 ++++++ src/Makefile | 2 +- src/ace.f90 | 15 +++- src/data_structures.f90 | 71 ++++++++++++++++- src/fileio.f90 | 141 ++++++++++++++++++++++++++++---- src/geometry.f90 | 54 ++++++++----- src/global.f90 | 38 ++++++++- src/input_sample | 9 ++- src/main.f90 | 21 ++--- src/output.f90 | 130 ++++++++++++++++++++++++++++++ src/physics.f90 | 173 +++++++++++++++++++++++++++++++++++++++- src/types.f90 | 42 +++++++--- 12 files changed, 645 insertions(+), 74 deletions(-) diff --git a/ChangeLog b/ChangeLog index 2e555e13a0..a9fd23cf08 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,26 @@ +2011-02-27 Paul Romano + + * ace.f90: Fixed minor bugs associated with reading angular + distributions. + * data_structures.f90: Added dict_keys routines. + * fileio.f90: Considerably changed read_count and read_cell to be + able to handle universes. + * geometry.f90: Added universe logic to find_cell and allowed it + to be recursively called. + * global.f90: Added universe array and dictionary, cell types, + pointer to lowest universe, angular distribution types, + interpolation types, as well as ERROR_CODE and adjusted str_to_int + to return ERROR_CODE on error. + * output.f90: Added print_cell, print_universe, and print_summary + subroutines. Need a few more useful routines like this! + * physics.f90: Added ability to sample cosine of scattering angle + from a 32-equiprobable bin distribution or a tabular + distribution. This is done by the sample_energy routine. The + rotate_angle subroutine takes a particles direction and applies + that scattering angle to it. Also added some logic for universes. + * types.f90: Changed some attributes in the Neutron structure, + added type to angular distribution within AceReaction. + 2011-02-22 Paul Romano * ace.f90: Added read_angular_dist and read_energy_dist diff --git a/src/Makefile b/src/Makefile index c8a737f0d8..359e34c36d 100644 --- a/src/Makefile +++ b/src/Makefile @@ -59,7 +59,7 @@ 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 -output.o: global.o +output.o: global.o types.o data_structures.o physics.o: types.o global.o mcnp_random.o geometry.o output.o \ search.o endf.o search.o: output.o diff --git a/src/ace.f90 b/src/ace.f90 index 53e55b03f3..bdd61a6779 100644 --- a/src/ace.f90 +++ b/src/ace.f90 @@ -365,6 +365,7 @@ contains NE = XSS(JXS9 + LOCB - 1) rxn % adist_n_energy = NE allocate(rxn % adist_energy(NE)) + allocate(rxn % adist_type(NE)) allocate(rxn % adist_location(NE)) ! read incoming energy grid and location of tables @@ -376,11 +377,16 @@ contains size = 0 do j = 1, NE LC = rxn % adist_location(j) - if (LC > 0) then + if (LC == 0) then + ! isotropic + rxn % adist_type(j) = ANGLE_ISOTROPIC + elseif (LC > 0) then ! 32 equiprobable bins - size = size + 32 + rxn % adist_type(j) = ANGLE_32_EQUI + size = size + 33 elseif (LC < 0) then ! tabular distribution + rxn % adist_type(j) = ANGLE_TABULAR NP = XSS(JXS9 + abs(LC)) size = size + 2 + 3*NP end if @@ -396,7 +402,10 @@ contains XSS_index = JXS9 + abs(LC) - 1 rxn % adist_data = get_real(size) -! print *, trim(table%name), rxn%MT, size + ! change location pointers since they are currently relative to + ! JXS(9) + LC = abs(rxn % adist_location(1)) + rxn % adist_location = abs(rxn % adist_location) - LC + 1 end do diff --git a/src/data_structures.f90 b/src/data_structures.f90 index 7462d793bb..01e6b19c7b 100644 --- a/src/data_structures.f90 +++ b/src/data_structures.f90 @@ -18,8 +18,8 @@ module data_structures implicit none - integer, parameter, private :: hash_size = 4993 - integer, parameter, private :: multiplier = 31 + integer, parameter :: hash_size = 4993 + integer, parameter :: multiplier = 31 integer, parameter :: DICT_NULL = 0 !===================================================================== @@ -82,6 +82,9 @@ module data_structures interface dict_hashkey module procedure dict_ci_hashkey, dict_ii_hashkey end interface + interface dict_keys + module procedure dict_ci_keys, dict_ii_keys + end interface contains @@ -1029,4 +1032,68 @@ contains end function dict_ii_hashkey +!===================================================================== +! DICT_CI_KEYS +!===================================================================== + + function dict_ci_keys(dict) result(head) + + type(DictionaryCI), pointer :: dict + type(ListKeyValueCI), pointer :: head + type(ListKeyValueCI), pointer :: current => null() + type(ListKeyValueCI), pointer :: elem => null() + + integer :: i, j + + head => null() + + do i = 1, size(dict%table) + elem => dict%table(i)%list + do while (associated(elem)) + if (.not. associated(head)) then + allocate(head) + current => head + else + allocate(current%next) + current => current%next + end if + current%data = elem%data + elem => elem%next + end do + end do + + end function dict_ci_keys + +!===================================================================== +! DICT_II_KEYS +!===================================================================== + + function dict_ii_keys(dict) result(head) + + type(DictionaryII), pointer :: dict + type(ListKeyValueII), pointer :: head + type(ListKeyValueII), pointer :: current => null() + type(ListKeyValueII), pointer :: elem => null() + + integer :: i, j + + head => null() + + do i = 1, size(dict%table) + elem => dict%table(i)%list + do while (associated(elem)) + if (.not. associated(head)) then + allocate(head) + current => head + else + allocate(current%next) + current => current%next + end if + current%data = elem%data + elem => elem%next + end do + end do + + end function dict_ii_keys + end module data_structures diff --git a/src/fileio.f90 b/src/fileio.f90 index 59f520bc1a..4c1cf81482 100644 --- a/src/fileio.f90 +++ b/src/fileio.f90 @@ -1,15 +1,22 @@ module fileio use global - use types, only: Cell, Surface, ExtSource + use types, only: Cell, Surface, ExtSource, ListKeyValueCI use string, only: split_string_wl, lower_case, str_to_real, split_string, & & concatenate use output, only: message, warning, error use data_structures, only: dict_create, dict_add_key, dict_get_key, & - & DICT_NULL + & dict_has_key, DICT_NULL, dict_keys, list_size implicit none + type(DictionaryII), pointer :: & ! this dictionary is used to count + & ucount_dict => null() ! the number of universes as well + ! as how many cells each universe + ! contains + + integer, allocatable :: index_cell_in_univ(:) + !===================================================================== ! READ_DATA interface allows data to be read with one function ! regardless of whether it is integer or real data. E.g. NXS and JXS @@ -60,8 +67,8 @@ contains !===================================================================== ! READ_COUNT makes a first pass through the input file to determine -! how many cells, surfaces, materials, sources, etc there are in the -! problem. +! how many cells, universes, surfaces, materials, sources, etc there +! are in the problem. !===================================================================== subroutine read_count(filename) @@ -73,14 +80,23 @@ contains integer :: in = 7 integer :: ioError integer :: n + integer :: count + integer :: universe_num + type(ListKeyValueII), pointer :: key_list + type(Universe), pointer :: univ + integer :: index msg = "First pass through input file..." call message(msg, 5) n_cells = 0 + n_universes = 0 n_surfaces = 0 n_materials = 0 + call dict_create(ucount_dict) + call dict_create(universe_dict) + open(FILE=filename, UNIT=in, STATUS='old', & & ACTION='read', IOSTAT=ioError) if (ioError /= 0) then @@ -96,6 +112,20 @@ contains select case (trim(words(1))) case ('cell') n_cells = n_cells + 1 + + ! For cells, we also need to check if there's a new universe + ! -- also for every cell add 1 to the count of cells for the + ! specified universe + universe_num = str_to_int(words(3)) + if (.not. dict_has_key(ucount_dict, universe_num)) then + n_universes = n_universes + 1 + count = 1 + call dict_add_key(universe_dict, universe_num, n_universes) + else + count = 1 + dict_get_key(ucount_dict, universe_num) + end if + call dict_add_key(ucount_dict, universe_num, count) + case ('surface') n_surfaces = n_surfaces + 1 case ('material') @@ -123,14 +153,47 @@ contains ! Allocate arrays for cells, surfaces, etc. allocate(cells(n_cells)) + allocate(universes(n_universes)) allocate(surfaces(n_surfaces)) allocate(materials(n_materials)) + ! Also allocate a list for keeping track of where cells have been + ! assigned in each universe + allocate(index_cell_in_univ(n_universes)) + index_cell_in_univ = 0 + ! Set up dictionaries call dict_create(cell_dict) call dict_create(surface_dict) call dict_create(material_dict) + ! We also need to allocate the cell count lists for each + ! universe. The logic for this is a little more convoluted. In + ! universe_dict, the (key,value) pairs are the uid of the universe + ! and the index in the array. In ucount_dict, it's the uid of the + ! universe and the number of cells. + + key_list => dict_keys(universe_dict) + do while (associated(key_list)) + ! find index of universe in universes array + index = key_list%data%value + univ => universes(index) + univ % uid = key_list%data%key + + ! check for lowest level universe + if (univ % uid == 0) BASE_UNIVERSE = index + + ! find cell count for this universe + count = dict_get_key(ucount_dict, key_list%data%key) + + ! allocate cell list for universe + allocate(univ % cells(count)) + univ % n_cells = count + + ! move to next universe + key_list => key_list%next + end do + end subroutine read_count !===================================================================== @@ -149,6 +212,7 @@ contains integer :: ioError integer :: n, i integer :: index_cell + integer :: index_universe integer :: index_surface integer :: index_material integer :: index_source @@ -157,6 +221,7 @@ contains call message(msg, 5) index_cell = 0 + index_universe = 0 index_surface = 0 index_material = 0 index_source = 0 @@ -205,6 +270,9 @@ contains end do + ! deallocate array + deallocate(index_cell_in_univ) + close(UNIT=in) call adjust_indices() @@ -245,13 +313,17 @@ contains end do ! adjust material indices - index = dict_get_key(material_dict, c%material) - if (index == DICT_NULL) then - msg = "Could not find material " // trim(int_to_str(c%material)) // & - & " specified on cell " // trim(int_to_str(c%uid)) - call error(msg) + if (c % fill /= 0) then + ! cell is filled with universe -- do nothing + else + index = dict_get_key(material_dict, c%material) + if (index == DICT_NULL) then + msg = "Could not find material " // trim(int_to_str(c%material)) // & + & " specified on cell " // trim(int_to_str(c%uid)) + call error(msg) + end if + c%material = index end if - c%material = index end do end subroutine adjust_indices @@ -268,32 +340,63 @@ contains integer :: ioError integer :: i + integer :: universe_num integer :: n_items character(250) :: msg character(32) :: word type(Cell), pointer :: this_cell => null() + type(Universe), pointer :: this_univ => null() this_cell => cells(index) ! Read cell identifier - read(words(2), FMT='(I8)', IOSTAT=ioError) this_cell%uid - if (ioError > 0) then + this_cell % uid = str_to_int(words(2)) + if (this_cell % uid == ERROR_CODE) then msg = "Invalid cell name: " // words(2) call error(msg) end if call dict_add_key(cell_dict, this_cell%uid, index) + + ! Read cell universe + universe_num = str_to_int(words(3)) + if (universe_num == ERROR_CODE) then + msg = "Invalid universe: " // words(3) + call error(msg) + end if + this_cell % universe = dict_get_key(universe_dict, universe_num) ! Read cell material - read(words(3), FMT='(I8)', IOSTAT=ioError) this_cell%material + if (trim(words(4)) == 'fill') then + this_cell % type = CELL_FILL + this_cell % material = 0 + n_items = n_words - 5 + + ! find universe + universe_num = str_to_int(words(5)) + if (universe_num == ERROR_CODE) then + msg = "Invalid universe fill: " // words(5) + call error(msg) + end if + this_cell % fill = dict_get_key(universe_dict, universe_num) + + else + this_cell % type = CELL_NORMAL + this_cell % material = str_to_int(words(4)) + this_cell % fill = 0 + if (this_cell % material == ERROR_CODE) then + msg = "Invalid material number: " // words(4) + call error(msg) + end if + n_items = n_words - 4 + end if ! Assign number of items - n_items = n_words - 3 this_cell%n_items = n_items ! Read list of surfaces allocate(this_cell%boundary_list(n_items)) do i = 1, n_items - word = words(i+3) + word = words(i+n_words-n_items) if (word(1:1) == '(') then this_cell%boundary_list(i) = OP_LEFT_PAREN elseif (word(1:1) == ')') then @@ -303,10 +406,16 @@ contains elseif (word(1:1) == '#') then this_cell%boundary_list(i) = OP_DIFFERENCE else - read(word, FMT='(I8)', IOSTAT=ioError) this_cell%boundary_list(i) + this_cell%boundary_list(i) = str_to_int(word) end if end do + ! Add cell to the cell list in the corresponding universe + i = this_cell % universe + this_univ => universes(i) + index_cell_in_univ(i) = index_cell_in_univ(i) + 1 + this_univ % cells(index_cell_in_univ(i)) = index + end subroutine read_cell !===================================================================== diff --git a/src/geometry.f90 b/src/geometry.f90 index 7fcc49e44a..1be471cd12 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -85,35 +85,47 @@ contains ! FIND_CELL determines what cell a source neutron is in !===================================================================== - subroutine find_cell(neut) + recursive subroutine find_cell(univ, neut, found) - type(Neutron), pointer :: neut + type(Universe), pointer :: univ + type(Neutron), pointer :: neut + logical, intent(inout) :: found - type(Cell), pointer :: this_cell - logical :: found_cell - character(250) :: msg - integer :: i + character(250) :: msg ! error message + integer :: i ! index over cells + type(Cell), pointer :: this_cell ! + type(Universe), pointer :: lower_univ + + found = .false. ! determine what region in - do i = 1, n_cells - this_cell => cells(i) + do i = 1, univ % n_cells + this_cell => cells(univ % cells(i)) + if (cell_contains(this_cell, neut)) then - neut%cell = i - found_cell = .true. + neut%cell = dict_get_key(cell_dict, this_cell % uid) - ! set current pointers - cCell => this_cell - cMaterial => materials(cCell%material) - exit + ! If this cell contains a universe of lattice, search for + ! the particle in that universe/lattice + if (this_cell % fill > 0) then + lower_univ => universes(this_cell % fill) + call find_cell(lower_univ, neut, found) + if (found) then + exit + else + msg = "Could not locate neutron in universe: " + call error(msg) + end if + else + ! set current pointers + found = .true. + cCell => this_cell + cMaterial => materials(cCell%material) + exit + end if end if - end do - ! if neutron couldn't be located, print error - if (.not. found_cell) then - write(msg, 100) "Could not locate cell for neutron at: ", neut%xyz -100 format (A,3ES11.3) - call error(msg) - end if + end do end subroutine find_cell diff --git a/src/global.f90 b/src/global.f90 index d645e172d2..0a9c15c93f 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -6,16 +6,19 @@ module global ! Main arrays for cells, surfaces, materials type(Cell), allocatable, target :: cells(:) + type(Universe), allocatable, target :: universes(:) type(Surface), allocatable, target :: surfaces(:) 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_surfaces ! # of surfaces integer :: n_materials ! # of materials ! These dictionaries provide a fast lookup mechanism type(DictionaryII), pointer :: cell_dict + type(DictionaryII), pointer :: universe_dict type(DictionaryII), pointer :: surface_dict type(DictionaryII), pointer :: material_dict type(DictionaryII), pointer :: isotope_dict @@ -30,6 +33,7 @@ module global ! Current cell, surface, material type(Cell), pointer :: cCell + type(Universe), pointer :: cUniverse type(Surface), pointer :: cSurface type(Material), pointer :: cMaterial @@ -79,6 +83,15 @@ module global & OP_UNION = huge(0) - 2, & ! Union operator & OP_DIFFERENCE = huge(0) - 3 ! Difference operator + ! Cell types + integer, parameter :: & + & CELL_NORMAL = 1, & ! Cell with a specified material + & CELL_FILL = 2, & ! Cell filled by a separate universe + & CELL_LATTICE = 3 ! Cell filled with a lattice + + ! array index of universe 0 + integer :: BASE_UNIVERSE + ! Surface types integer, parameter :: & & SURF_PX = 1, & ! Plane parallel to x-plane @@ -112,12 +125,28 @@ module global & PROB_SOURCE = 1, & ! External source problem & PROB_CRITICALITY = 2 ! Criticality problem + ! Angular distribution type + integer, parameter :: & + & ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution + & ANGLE_32_EQUI = 2, & ! 32 equiprobable bins + & ANGLE_TABULAR = 3 ! Tabular angular distribution + + ! Interpolation flag + integer, parameter :: & + & HISTOGRAM = 1, & ! Histogram interpolation on angle + & LINEARLINEAR = 2 ! Linear-linear interpolation on angle + ! Particle type integer, parameter :: & & NEUTRON_ = 1, & & PHOTON_ = 2, & & ELECTRON_ = 3 + ! Integer code for read error -- better hope this number is never + ! used in an input file! + integer, parameter :: & + & ERROR_CODE = -huge(0) + integer :: verbosity = 5 integer, parameter :: max_words = 500 integer, parameter :: max_line = 250 @@ -192,15 +221,18 @@ contains integer :: num character(5) :: fmt - integer :: w + integer :: w + integer :: ioError + ! Determine width of string w = len_trim(str) ! Create format specifier for reading string - write(fmt, '("(I",I2,")")') w + write(UNIT=fmt, FMT='("(I",I2,")")') w ! read string into integer - read(str,fmt) num + read(UNIT=str, FMT=fmt, IOSTAT=ioError) num + if (ioError > 0) num = ERROR_CODE end function str_to_int diff --git a/src/input_sample b/src/input_sample index be9c70d37d..e067a1dd84 100644 --- a/src/input_sample +++ b/src/input_sample @@ -1,10 +1,13 @@ # test input file -cell 100 40 -1 -cell 200 41 1 -2 +cell 1 0 fill 37 -2 +cell 100 37 40 -1 +cell 101 37 41 1 +cell 2 0 40 2 -3 surface 1 sph 0 0 0 3 surface 2 sph 0 0 0 5 +surface 3 sph 0 0 0 8 material 40 -4.5 & 92238.03c 1.0 @@ -13,7 +16,7 @@ material 41 -1.0 & 1001.03c 2.0 & 8016.03c 1.0 -source box -1 -1 -1 1 1 1 +source box -4 -4 -4 4 4 4 xs_library endfb7 xs_data /opt/serpent/xsdata/endfb7 diff --git a/src/main.f90 b/src/main.f90 index df74c88bb6..771828ae5a 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -1,17 +1,17 @@ program main use global - use fileio, only: read_input, read_command_line, read_count, & - & normalize_ao - use output, only: title, echo_input, message, warning, error - 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 - use physics, only: transport - use data_structures, only: dict_create, dict_add_key, dict_get_key + use fileio, only: read_input, read_command_line, read_count, & + & normalize_ao + use output, only: title, echo_input, message, warning, error, & + & print_summary + 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 + use physics, only: transport use cross_section, only: read_xsdata, material_total_xs - use ace, only: read_xs - use energy_grid, only: unionized_grid, original_indices + use ace, only: read_xs + use energy_grid, only: unionized_grid, original_indices implicit none @@ -60,6 +60,7 @@ program main call material_total_xs() call echo_input() + call print_summary() ! create source particles call init_source() diff --git a/src/output.f90 b/src/output.f90 index e864ca432b..ad74f54803 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -2,6 +2,8 @@ module output use ISO_FORTRAN_ENV use global + use types, only: Cell, Universe, Surface + use data_structures, only: dict_get_key implicit none @@ -218,4 +220,132 @@ module output end subroutine get_today +!===================================================================== +! PRINT_CELL displays the attributes of a cell +!===================================================================== + + subroutine print_cell(c) + + type(Cell), pointer :: c + + integer :: ou + integer :: temp + integer :: i + type(Universe), pointer :: u => null() + type(Material), pointer :: m => null() + character(250) :: string + + ou = OUTPUT_UNIT + + write(ou,*) 'Cell ' // int_to_str(c % uid) + temp = dict_get_key(cell_dict, c % uid) + write(ou,*) ' Array Index = ' // int_to_str(temp) + u => universes(c % universe) + write(ou,*) ' Universe = ' // int_to_str(u % uid) + if (c % fill == 0) then + write(ou,*) ' Fill = NONE' + else + u => universes(c % fill) + write(ou,*) ' Fill = ' // int_to_str(u % uid) + end if + if (c % material == 0) then + write(ou,*) ' Material = NONE' + else + m => materials(c % material) + write(ou,*) ' Material = ' // int_to_str(m % uid) + end if + string = "" + do i = 1, c % n_items + select case (c % boundary_List(i)) + case (OP_LEFT_PAREN) + string = trim(string) // ' (' + case (OP_RIGHT_PAREN) + string = trim(string) // ' )' + case (OP_UNION) + string = trim(string) // ' :' + case (OP_DIFFERENCE) + string = trim(string) // ' !' + case default + string = trim(string) // ' ' // int_to_str(c % boundary_list(i)) + end select + end do + write(ou,*) ' Surface Specification:' // trim(string) + write(ou,*) + + ! nullify associated pointers + nullify(u) + nullify(m) + + end subroutine print_cell + +!===================================================================== +! PRINT_UNIVERSE displays the attributes of a universe +!===================================================================== + + subroutine print_universe(univ) + + type(Universe), pointer :: univ + + integer :: ou + integer :: i + character(250) :: string + type(Cell), pointer :: c + + ou = OUTPUT_UNIT + + write(ou,*) 'Universe ' // int_to_str(univ % uid) + write(ou,*) ' Level = ' // int_to_str(univ % level) + string = "" + do i = 1, univ % n_cells + c => cells(univ % cells(i)) + string = trim(string) // ' ' // int_to_str(c % uid) + end do + write(ou,*) ' Cells =' // trim(string) + write(ou,*) + + end subroutine print_universe + +!===================================================================== +! PRINT_SUMMARY displays the attributes of all cells, universes, +! surfaces and materials read in the input file. Very useful for +! debugging! +!===================================================================== + + subroutine print_summary() + + type(Surface), pointer :: s + type(Cell), pointer :: c + type(Universe), pointer :: u + type(Material), pointer :: m + integer :: i + integer :: ou + + ou = OUTPUT_UNIT + + ! print summary of cells + write(ou,*) '=============================================' + write(ou,*) '=> CELL SUMMARY <=' + write(ou,*) '=============================================' + write(ou,*) + do i = 1, n_cells + c => cells(i) + call print_cell(c) + end do + + ! print summary of universes + write(ou,*) '=============================================' + write(ou,*) '=> UNIVERSE SUMMARY <=' + write(ou,*) '=============================================' + write(ou,*) + do i = 1, n_universes + u => universes(i) + call print_universe(u) + end do + + ! print summary of surfaces + + ! print summary of materials + + end subroutine print_summary + end module output diff --git a/src/physics.f90 b/src/physics.f90 index 99f1bb209d..6786dd982f 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -4,7 +4,7 @@ module physics use geometry, only: find_cell, dist_to_boundary, cross_boundary use types, only: Neutron use mcnp_random, only: rang - use output, only: error, message + use output, only: error, message, print_universe use search, only: binary_search use endf, only: reaction_name @@ -19,11 +19,12 @@ contains subroutine transport(neut) - type(Neutron), pointer, intent(inout) :: neut + type(Neutron), pointer :: neut character(250) :: msg integer :: i, surf logical :: alive + logical :: found_cell real(8) :: d_to_boundary real(8) :: d_to_collision @@ -32,11 +33,23 @@ contains real(8) :: f ! interpolation factor real(8) :: tmp(3) integer :: IE ! index on energy grid + type(Universe), pointer :: univ - ! determine what cell the particle is in if (neut%cell == 0) then - call find_cell(neut) + univ => universes(BASE_UNIVERSE) + call find_cell(univ, neut, found_cell) + print *, sqrt(neut%xyz(1)**2 + neut%xyz(2)**2 + neut%xyz(3)**2), cells(neut%cell) % uid + + ! if neutron couldn't be located, print error + if (.not. found_cell) then + write(msg, 100) "Could not locate cell for neutron at: ", neut%xyz +100 format (A,3ES11.3) + call error(msg) + end if end if + + return + if (verbosity >= 10) then msg = "=== Particle " // trim(int_to_str(neut%uid)) // " ===" call message(msg, 10) @@ -268,6 +281,158 @@ contains end subroutine n_gamma +!===================================================================== +! SAMPLE_ANGLE samples the cosine of the angle between incident and +! exiting particle directions either from 32 equiprobable bins or from +! a tabular distribution. +!===================================================================== + + subroutine sample_angle(rxn, E) + + type(AceReaction), pointer :: rxn ! reaction + real(8), intent(in) :: E ! incoming energy + + real(8) :: xi ! random number on [0,1) + integer :: interp ! type of interpolation + integer :: type ! angular distribution type + integer :: i ! incoming energy bin + integer :: n ! number of incoming energy bins + integer :: loc ! location in data array + integer :: np ! number of points in cos distribution + integer :: bin ! cosine bin + real(8) :: f ! interpolation factor + real(8) :: mu0 ! cosine in bin b + real(8) :: mu1 ! cosine in bin b+1 + real(8) :: mu ! final cosine sampled + real(8) :: c ! cumulative distribution frequency + real(8) :: p0,p1 ! probability distribution + character(250) :: msg ! error message + + ! determine number of incoming energies + n = rxn % adist_n_energy + + ! find energy bin and calculate interpolation factor -- if the + ! energy is outside the range of the tabulated energies, choose + ! the first or last bins + if (E < rxn%adist_energy(1)) then + i = 1 + f = 0.0 + elseif (E > rxn%adist_energy(n)) then + i = n - 1 + f = 1.0 + else + i = binary_search(rxn%adist_energy, n, E) + f = (E - rxn % adist_energy(i)) / & + & (rxn % adist_energy(i+1) - rxn % adist_energy(i)) + end if + + ! Sample between the ith and (i+1)th bin + if (f > rang()) i = i + 1 + + ! check whether this is a 32-equiprobable bin or a tabular + ! distribution + loc = rxn % adist_location(i) + type = rxn % adist_type(i) + if (type == ANGLE_ISOTROPIC) then + mu = 2.0_8 * rang() - 1 + elseif (type == ANGLE_32_EQUI) then + ! sample cosine bin + xi = rang() + bin = 1 + int(32.0_8*xi) + + ! calculate cosine + mu0 = rxn % adist_data(loc + bin - 1) + mu1 = rxn % adist_data(loc + bin) + mu = mu0 + (32.0_8 * xi - bin) * (mu1 - mu0) + + elseif (type == ANGLE_TABULAR) then + interp = rxn % adist_data(loc) + np = rxn % adist_data(loc+1) + + ! determine outgoing cosine bin + xi = rang() + do bin = loc+2, loc+1+NP + c = rxn % adist_data(bin+2*np) + if (xi > c) exit + end do + + p0 = rxn % adist_data(bin+np) + mu0 = rxn % adist_data(bin) + if (interp == HISTOGRAM) then + ! Histogram interpolation + mu = mu0 + (xi - c)/p0 + + elseif (interp == LINEARLINEAR) then + ! Linear-linear interpolation -- not sure how you come about + ! the formula given in the MCNP manual + p1 = rxn % adist_data(bin+np+1) + mu1 = rxn % adist_data(bin+1) + + f = (p1 - p0)/(mu1 - mu0) + if (f == ZERO) then + mu = mu0 + (xi - c)/p0 + else + mu = mu0 + (sqrt(p0*p0 + 2*f*(xi - c))-p0)/f + end if + else + msg = "Unknown interpolation type: " // trim(int_to_str(interp)) + call error(msg) + end if + + else + msg = "Unknown angular distribution type: " // trim(int_to_str(type)) + call error(msg) + end if + + end subroutine sample_angle + +!===================================================================== +! ROTATE_ANGLE rotates direction cosines through a polar angle whose +! cosine is mu and through an azimuthal angle sampled uniformly. Note +! that this is done with direct sampling rather than rejection as is +! done in MCNP and SERPENT. +!===================================================================== + + subroutine rotate_angle(mu, u, v, w) + + real(8), intent(in) :: mu ! cosine of angle + real(8), intent(inout) :: u + real(8), intent(inout) :: v + real(8), intent(inout) :: w + + real(8) :: phi, sinphi, cosphi + real(8) :: a,b + real(8) :: u0, v0, w0 + + ! Copy original directional cosines + u0 = u + v0 = v + w0 = w + + ! Sample azimuthal angle in [0,2pi) + phi = 2.0_8 * PI * rang() + + ! Precompute factors to save flops + sinphi = sin(phi) + cosphi = cos(phi) + a = sqrt(1.0_8 - mu*mu) + b = sqrt(1.0_8 - w*w) + + ! Need to treat special case where sqrt(1 - w**2) is close to zero + ! by expanding about the v component rather than the w component + if (b > 1e-10) then + u = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b + v = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b + w = mu*w0 - a*b*cosphi + else + b = sqrt(1.0_8 - v*v) + u = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b + v = mu*v0 - a*b*cosphi + w = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b + end if + + end subroutine rotate_angle + !===================================================================== ! MAXWELL_SPECTRUM samples an energy from the Maxwell fission ! distribution based on a rejection sampling scheme. This is described diff --git a/src/types.f90 b/src/types.f90 index b5ed76f187..ceaa2ad6a2 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -2,6 +2,21 @@ module types implicit none +!===================================================================== +! UNIVERSE +!===================================================================== + + type Universe + integer :: uid + integer :: type + integer :: level + integer :: n_cells + integer, allocatable :: cells(:) + real(8) :: x0 + real(8) :: y0 + real(8) :: z0 + end type Universe + !===================================================================== ! SURFACE type defines a first- or second-order surface that can be ! used to construct closed volumes (cells) @@ -22,9 +37,12 @@ module types type Cell integer :: uid + integer :: type + integer :: universe + integer :: fill + integer :: material integer :: n_items integer, allocatable :: boundary_list(:) - integer :: material end type Cell !===================================================================== @@ -33,16 +51,17 @@ module types !===================================================================== type Neutron - integer :: uid ! Unique ID - real(8) :: xyz(3) ! location - real(8) :: uvw(3) ! directional cosines - real(8) :: E ! energy - integer :: IE ! index on energy grid - real(8) :: interp ! interpolation factor for energy grid - integer :: cell ! current cell - integer :: surface ! current surface - real(8) :: wgt ! particle weight - logical :: alive ! is particle alive? + integer :: uid ! Unique ID + real(8) :: xyz(3) ! location + real(8) :: uvw(3) ! directional cosines + real(8) :: E ! energy + integer :: IE ! index on energy grid + real(8) :: interp ! interpolation factor for energy grid + integer :: cell ! current cell + integer :: universe ! current universe + integer :: surface ! current surface + real(8) :: wgt ! particle weight + logical :: alive ! is particle alive? end type Neutron !===================================================================== @@ -115,6 +134,7 @@ module types ! Secondary angle distribution integer :: adist_n_energy ! # of incoming energies real(8), allocatable :: adist_energy(:) ! incoming energy grid + integer, allocatable :: adist_type(:) ! type of distribution integer, allocatable :: adist_location(:) ! location of each table real(8), allocatable :: adist_data(:) ! angular distribution data