From 1c42639ef5da861c7cf7216964bcf10994a4c6df Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 10 Mar 2011 22:37:17 +0000 Subject: [PATCH] Changed attribute names in Cell, Neutron is now Particle. --- ChangeLog | 17 ++++ src/fileio.f90 | 64 ++++++------ src/geometry.f90 | 260 ++++++++++++++++++++++++----------------------- src/global.f90 | 10 +- src/main.f90 | 8 +- src/output.f90 | 6 +- src/physics.f90 | 116 ++++++++++----------- src/source.f90 | 4 +- src/types.f90 | 63 ++++++------ 9 files changed, 286 insertions(+), 262 deletions(-) diff --git a/ChangeLog b/ChangeLog index 41f0a22ec3..7440016236 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,20 @@ +2011-03-10 Paul Romano + + * types.f90: The major changes are that the n_items and + boundary_list attributes of 'Cell' are now n_surfaces and + surfaces. Even though this doesn't quite make sense since + parentheses and union operators are not surfaces per-se, it is + easier to understand within other modules with these names. The + 'Neutron' type is now called 'Particle' and will be used for all + different particles (eventually). + * fileio.f90: Changed aforementioned names. + * geometry.f90: Changed aforementioned names. + * global.f90: Changed aforementioned names. + * main.f90: Changed aforementioned names. + * output.f90: Changed aforementioned names. + * physics.f90: Changed aforementioned names. + * source.f90: Changed aforementioned names. + 2011-03-04 Paul Romano * fileio.f90: Added build_universe routine which sets pointers diff --git a/src/fileio.f90 b/src/fileio.f90 index 7430af6c89..73dfbe0700 100644 --- a/src/fileio.f90 +++ b/src/fileio.f90 @@ -305,8 +305,8 @@ contains do i = 1, n_cells ! adjust boundary list c => cells(i) - do j = 1, c%n_items - surf_num = c%boundary_list(j) + do j = 1, c%n_surfaces + surf_num = c%surfaces(j) if (surf_num < OP_DIFFERENCE) then index = dict_get_key(surface_dict, abs(surf_num)) if (index == DICT_NULL) then @@ -315,7 +315,7 @@ contains & " specified on cell " // trim(int_to_str(c%uid)) call error(msg) end if - c%boundary_list(j) = sign(index, surf_num) + c%surfaces(j) = sign(index, surf_num) end if end do @@ -383,21 +383,20 @@ contains integer :: ioError integer :: i integer :: universe_num - integer :: n_items + integer :: n_surfaces character(250) :: msg - character(32) :: word - type(Cell), pointer :: this_cell => null() - type(Universe), pointer :: this_univ => null() + character(32) :: word + type(Cell), pointer :: c => null() - this_cell => cells(index) + c => cells(index) ! Read cell identifier - this_cell % uid = str_to_int(words(2)) - if (this_cell % uid == ERROR_CODE) then + c % uid = str_to_int(words(2)) + if (c % 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) + call dict_add_key(cell_dict, c%uid, index) ! Read cell universe universe_num = str_to_int(words(3)) @@ -405,13 +404,13 @@ contains msg = "Invalid universe: " // words(3) call error(msg) end if - this_cell % universe = dict_get_key(universe_dict, universe_num) + c % universe = dict_get_key(universe_dict, universe_num) ! Read cell material if (trim(words(4)) == 'fill') then - this_cell % type = CELL_FILL - this_cell % material = 0 - n_items = n_words - 5 + c % type = CELL_FILL + c % material = 0 + n_surfaces = n_words - 5 ! find universe universe_num = str_to_int(words(5)) @@ -419,44 +418,43 @@ contains msg = "Invalid universe fill: " // words(5) call error(msg) end if - this_cell % fill = dict_get_key(universe_dict, universe_num) + c % 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 + c % type = CELL_NORMAL + c % material = str_to_int(words(4)) + c % fill = 0 + if (c % material == ERROR_CODE) then msg = "Invalid material number: " // words(4) call error(msg) end if - n_items = n_words - 4 + n_surfaces = n_words - 4 end if ! Assign number of items - this_cell%n_items = n_items + c%n_surfaces = n_surfaces ! Read list of surfaces - allocate(this_cell%boundary_list(n_items)) - do i = 1, n_items - word = words(i+n_words-n_items) + allocate(c%surfaces(n_surfaces)) + do i = 1, n_surfaces + word = words(i+n_words-n_surfaces) if (word(1:1) == '(') then - this_cell%boundary_list(i) = OP_LEFT_PAREN + c % surfaces(i) = OP_LEFT_PAREN elseif (word(1:1) == ')') then - this_cell%boundary_list(i) = OP_RIGHT_PAREN + c % surfaces(i) = OP_RIGHT_PAREN elseif (word(1:1) == ':') then - this_cell%boundary_list(i) = OP_UNION + c % surfaces(i) = OP_UNION elseif (word(1:1) == '#') then - this_cell%boundary_list(i) = OP_DIFFERENCE + c % surfaces(i) = OP_DIFFERENCE else - this_cell%boundary_list(i) = str_to_int(word) + c % surfaces(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) + i = c % universe index_cell_in_univ(i) = index_cell_in_univ(i) + 1 - this_univ % cells(index_cell_in_univ(i)) = index + universes(i) % cells(index_cell_in_univ(i)) = index end subroutine read_cell diff --git a/src/geometry.f90 b/src/geometry.f90 index 0bb71b0928..483d61c042 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -13,28 +13,28 @@ contains ! CELL_CONTAINS determines whether a given point is inside a cell !===================================================================== - function cell_contains(c, neut) result(in_cell) + function cell_contains(c, p) result(in_cell) - type(Cell), pointer :: c - type(Neutron), pointer :: neut - logical :: in_cell + type(Cell), pointer :: c + type(Particle), pointer :: p + logical :: in_cell integer, allocatable :: expression(:) integer :: specified_sense integer :: actual_sense - integer :: n_boundaries + integer :: n_surfaces integer :: i, j integer :: surf_num integer :: current_surface type(Surface), pointer :: surf => null() character(250) :: msg - current_surface = neut%surface + current_surface = p%surface - n_boundaries = size(c%boundary_list) - allocate(expression(n_boundaries)) - expression = c%boundary_list - do i = 1, n_boundaries + n_surfaces = size(c%surfaces) + allocate(expression(n_surfaces)) + expression = c%surfaces + do i = 1, n_surfaces ! Don't change logical operator if (expression(i) >= OP_DIFFERENCE) then @@ -47,11 +47,11 @@ contains ! Check if the particle is currently on the specified surface if (surf_num == current_surface) then - ! neutron is on specified surface heading into cell + ! particle is on specified surface heading into cell expression(i) = 1 cycle elseif (surf_num == -current_surface) then - ! neutron is on specified surface, but heading other + ! particle is on specified surface, but heading other ! direction expression(i) = 0 cycle @@ -59,7 +59,7 @@ contains ! Compare sense of point to specified sense specified_sense = sign(1,expression(i)) - call sense(surf, neut%xyz, actual_sense) + call sense(surf, p%xyz, actual_sense) if (actual_sense == specified_sense) then expression(i) = 1 else @@ -82,44 +82,48 @@ contains end function cell_contains !===================================================================== -! FIND_CELL determines what cell a source neutron is in +! FIND_CELL determines what cell a source particle is in within a +! particular universe. If the base universe is passed, the particle +! should be found as long as it's within the geometry !===================================================================== - recursive subroutine find_cell(univ, neut, found) + recursive subroutine find_cell(univ, p, found) - type(Universe), pointer :: univ - type(Neutron), pointer :: neut - logical, intent(inout) :: found + type(Universe), pointer :: univ ! universe to search in + type(Particle), pointer :: p ! pointer to particle + logical, intent(inout) :: found ! particle found? character(250) :: msg ! error message integer :: i ! index over cells - type(Cell), pointer :: this_cell ! - type(Universe), pointer :: lower_univ + type(Cell), pointer :: c ! pointer to cell + type(Universe), pointer :: lower_univ ! if particle is in lower + ! universe, use this pointer + ! to call recursively found = .false. ! determine what region in do i = 1, univ % n_cells - this_cell => cells(univ % cells(i)) + c => cells(univ % cells(i)) - if (cell_contains(this_cell, neut)) then - neut%cell = dict_get_key(cell_dict, this_cell % uid) + if (cell_contains(c, p)) then + p%cell = dict_get_key(cell_dict, c % uid) ! 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 (c % fill > 0) then + lower_univ => universes(c % fill) + call find_cell(lower_univ, p, found) if (found) then exit else - msg = "Could not locate neutron in universe: " + msg = "Could not locate particle in universe: " call error(msg) end if else ! set current pointers found = .true. - cCell => this_cell + cCell => c cMaterial => materials(cCell%material) exit end if @@ -130,12 +134,12 @@ contains end subroutine find_cell !===================================================================== -! CROSS_BOUNDARY moves a neutron into a new cell +! CROSS_BOUNDARY moves a particle into a new cell !===================================================================== - subroutine cross_boundary(neut) + subroutine cross_boundary(p) - type(Neutron), pointer :: neut + type(Particle), pointer :: p type(Surface), pointer :: surf type(Cell), pointer :: c @@ -145,7 +149,7 @@ contains logical :: found type(Universe), pointer :: lower_univ => null() - surf => surfaces(abs(neut%surface)) + surf => surfaces(abs(p%surface)) if (verbosity >= 10) then msg = " Crossing surface " // trim(int_to_str(surf%uid)) call message(msg, 10) @@ -153,7 +157,7 @@ contains ! check for leakage if (surf%bc == BC_VACUUM) then - neut%alive = .false. + p%alive = .false. if (verbosity >= 10) then msg = " Leaked out of surface " // trim(int_to_str(surf%uid)) call message(msg, 10) @@ -161,46 +165,46 @@ contains return end if - if (neut%surface > 0 .and. allocated(surf%neighbor_pos)) then + if (p%surface > 0 .and. allocated(surf%neighbor_pos)) then ! If coming from negative side of surface, search all the ! neighboring cells on the positive side do i = 1, size(surf%neighbor_pos) index_cell = surf%neighbor_pos(i) c => cells(index_cell) - if (cell_contains(c, neut)) then + if (cell_contains(c, p)) then if (c % fill > 0) then lower_univ => universes(c % fill) - call find_cell(lower_univ, neut, found) + call find_cell(lower_univ, p, found) if (.not. found) then - msg = "Could not locate neutron in universe: " + msg = "Could not locate particle in universe: " call error(msg) end if else ! set current pointers - neut%cell = index_cell + p%cell = index_cell cCell => c cMaterial => materials(cCell%material) end if return end if end do - elseif (neut%surface < 0 .and. allocated(surf%neighbor_neg)) then + elseif (p%surface < 0 .and. allocated(surf%neighbor_neg)) then ! If coming from positive side of surface, search all the ! neighboring cells on the negative side do i = 1, size(surf%neighbor_neg) index_cell = surf%neighbor_neg(i) c => cells(index_cell) - if (cell_contains(c, neut)) then + if (cell_contains(c, p)) then if (c % fill > 0) then lower_univ => universes(c % fill) - call find_cell(lower_univ, neut, found) + call find_cell(lower_univ, p, found) if (.not. found) then - msg = "Could not locate neutron in universe: " + msg = "Could not locate particle in universe: " call error(msg) end if else ! set current pointers - neut%cell = index_cell + p%cell = index_cell cCell => c cMaterial => materials(cCell%material) end if @@ -213,8 +217,8 @@ contains ! cells do i = 1, size(cells) c => cells(i) - if (cell_contains(c, neut)) then - neut%cell = i + if (cell_contains(c, p)) then + p%cell = i cCell => c cMaterial => materials(cCell%material) return @@ -222,7 +226,7 @@ contains end do ! Couldn't find next cell anywhere! - msg = "After neutron crossed surface " // trim(int_to_str(neut%surface)) // & + msg = "After particle crossed surface " // trim(int_to_str(p%surface)) // & & ", it could not be located in any cell and it did not leak." call error(msg) @@ -230,14 +234,14 @@ contains !===================================================================== ! DIST_TO_BOUNDARY calculates the distance to the nearest boundary of -! the cell 'cl' for a particle 'neut' traveling in a certain +! the cell 'cl' for a particle 'p' traveling in a certain ! direction. For a cell in a subuniverse that has a parent cell, also ! include the surfaces of the edge of the universe. !===================================================================== - subroutine dist_to_boundary(neut, dist, surf, other_cell) + subroutine dist_to_boundary(p, dist, surf, other_cell) - type(Neutron), pointer :: neut + type(Particle), pointer :: p real(8), intent(out) :: dist integer, intent(out) :: surf integer, optional, intent(in) :: other_cell @@ -246,7 +250,7 @@ contains type(Cell), pointer :: parent_p => null() type(Surface), pointer :: surf_p => null() integer :: i - integer :: n_boundaries, n1, n2 + integer :: n_surfaces, n1, n2 integer, allocatable :: expression(:) integer :: index_surf integer :: current_surf @@ -261,36 +265,36 @@ contains if (present(other_cell)) then cell_p => cells(other_cell) else - cell_p => cells(neut%cell) + cell_p => cells(p%cell) end if - current_surf = neut%surface + current_surf = p%surface ! determine number of surfaces to check - n1 = cell_p % n_items + n1 = cell_p % n_surfaces n2 = 0 if (cell_p % parent > 0) then parent_p => cells(cell_p % parent) - n2 = parent_p % n_items + n2 = parent_p % n_surfaces end if - n_boundaries = n1 + n2 + n_surfaces = n1 + n2 ! allocate space and assign expression - allocate(expression(n_boundaries)) - expression(1:n1) = cell_p%boundary_list + allocate(expression(n_surfaces)) + expression(1:n1) = cell_p % surfaces if (cell_p % parent > 0) then - expression(n1+1:n1+n2) = parent_p % boundary_list + expression(n1+1:n1+n2) = parent_p % surfaces end if ! loop over all surfaces dist = INFINITY - do i = 1, n_boundaries - x = neut%xyz(1) - y = neut%xyz(2) - z = neut%xyz(3) - u = neut%uvw(1) - v = neut%uvw(2) - w = neut%uvw(3) + do i = 1, n_surfaces + x = p % xyz(1) + y = p % xyz(2) + z = p % xyz(3) + u = p % uvw(1) + v = p % uvw(2) + w = p % uvw(3) ! check for coincident surface index_surf = expression(i) @@ -306,7 +310,7 @@ contains if (u == ZERO) then d = INFINITY else - x0 = surf_p%coeffs(1) + x0 = surf_p % coeffs(1) d = (x0 - x)/u if (d < ZERO) d = INFINITY end if @@ -315,7 +319,7 @@ contains if (v == ZERO) then d = INFINITY else - y0 = surf_p%coeffs(1) + y0 = surf_p % coeffs(1) d = (y0 - y)/v if (d < ZERO) d = INFINITY end if @@ -324,16 +328,16 @@ contains if (w == ZERO) then d = INFINITY else - z0 = surf_p%coeffs(1) + z0 = surf_p % coeffs(1) d = (z0 - z)/w if (d < ZERO) d = INFINITY end if case (SURF_PLANE) - A = surf_p%coeffs(1) - B = surf_p%coeffs(2) - C = surf_p%coeffs(3) - D = surf_p%coeffs(4) + A = surf_p % coeffs(1) + B = surf_p % coeffs(2) + C = surf_p % coeffs(3) + D = surf_p % coeffs(4) tmp = A*u + B*v + C*w if (tmp == ZERO) then @@ -348,9 +352,9 @@ contains if (a == ZERO) then d = INFINITY else - y0 = surf_p%coeffs(1) - z0 = surf_p%coeffs(2) - r = surf_p%coeffs(3) + y0 = surf_p % coeffs(1) + z0 = surf_p % coeffs(2) + r = surf_p % coeffs(3) y = y - y0 z = z - z0 @@ -383,9 +387,9 @@ contains if (a == ZERO) then d = INFINITY else - x0 = surf_p%coeffs(1) - z0 = surf_p%coeffs(2) - r = surf_p%coeffs(3) + x0 = surf_p % coeffs(1) + z0 = surf_p % coeffs(2) + r = surf_p % coeffs(3) x = x - x0 z = z - z0 @@ -418,9 +422,9 @@ contains if (a == ZERO) then d = INFINITY else - x0 = surf_p%coeffs(1) - y0 = surf_p%coeffs(2) - r = surf_p%coeffs(3) + x0 = surf_p % coeffs(1) + y0 = surf_p % coeffs(2) + r = surf_p % coeffs(3) x = x - x0 y = y - y0 @@ -454,10 +458,10 @@ contains end if case (SURF_SPHERE) - x0 = surf_p%coeffs(1) - y0 = surf_p%coeffs(2) - z0 = surf_p%coeffs(3) - r = surf_p%coeffs(4) + x0 = surf_p % coeffs(1) + y0 = surf_p % coeffs(2) + z0 = surf_p % coeffs(3) + r = surf_p % coeffs(4) x = x - x0 y = y - y0 @@ -531,56 +535,56 @@ contains y = xyz(2) z = xyz(3) - select case (surf%type) + select case (surf % type) case (SURF_PX) - x0 = surf%coeffs(1) + x0 = surf % coeffs(1) func = x - x0 case (SURF_PY) - y0 = surf%coeffs(1) + y0 = surf % coeffs(1) func = y - y0 case (SURF_PZ) - z0 = surf%coeffs(1) + z0 = surf % coeffs(1) func = z - z0 case (SURF_PLANE) - A = surf%coeffs(1) - B = surf%coeffs(2) - C = surf%coeffs(3) - D = surf%coeffs(4) + A = surf % coeffs(1) + B = surf % coeffs(2) + C = surf % coeffs(3) + D = surf % coeffs(4) func = A*x + B*y + C*z - D case (SURF_CYL_X) - y0 = surf%coeffs(1) - z0 = surf%coeffs(2) - r = surf%coeffs(3) + y0 = surf % coeffs(1) + z0 = surf % coeffs(2) + r = surf % coeffs(3) func = (y-y0)**2 + (z-z0)**2 - r**2 case (SURF_CYL_Y) - x0 = surf%coeffs(1) - z0 = surf%coeffs(2) - r = surf%coeffs(3) + x0 = surf % coeffs(1) + z0 = surf % coeffs(2) + r = surf % coeffs(3) func = (x-x0)**2 + (z-z0)**2 - r**2 case (SURF_CYL_Z) - x0 = surf%coeffs(1) - y0 = surf%coeffs(2) - r = surf%coeffs(3) + x0 = surf % coeffs(1) + y0 = surf % coeffs(2) + r = surf % coeffs(3) func = (x-x0)**2 + (y-y0)**2 - r**2 case (SURF_SPHERE) - x0 = surf%coeffs(1) - y0 = surf%coeffs(2) - z0 = surf%coeffs(3) - r = surf%coeffs(4) + x0 = surf % coeffs(1) + y0 = surf % coeffs(2) + z0 = surf % coeffs(3) + r = surf % coeffs(4) func = (x-x0)**2 + (y-y0)**2 + (z-z0)**2 - r**2 case (SURF_BOX_X) - y0 = surf%coeffs(1) - z0 = surf%coeffs(2) - y1 = surf%coeffs(3) - z1 = surf%coeffs(4) + y0 = surf % coeffs(1) + z0 = surf % coeffs(2) + y1 = surf % coeffs(3) + z1 = surf % coeffs(4) if (y >= y0 .and. y < y1 .and. z >= z0 .and. z < z1) then s = SENSE_NEGATIVE else @@ -589,10 +593,10 @@ contains return case (SURF_BOX_Y) - x0 = surf%coeffs(1) - z0 = surf%coeffs(2) - x1 = surf%coeffs(3) - z1 = surf%coeffs(4) + x0 = surf % coeffs(1) + z0 = surf % coeffs(2) + x1 = surf % coeffs(3) + z1 = surf % coeffs(4) if (x >= x0 .and. x < x1 .and. z >= z0 .and. z < z1) then s = SENSE_NEGATIVE else @@ -601,10 +605,10 @@ contains return case (SURF_BOX_Z) - x0 = surf%coeffs(1) - y0 = surf%coeffs(2) - x1 = surf%coeffs(3) - y1 = surf%coeffs(4) + x0 = surf % coeffs(1) + y0 = surf % coeffs(2) + x1 = surf % coeffs(3) + y1 = surf % coeffs(4) if (x >= x0 .and. x < x1 .and. y >= y0 .and. y < y1) then s = SENSE_NEGATIVE else @@ -613,12 +617,12 @@ contains return case (SURF_BOX) - x0 = surf%coeffs(1) - y0 = surf%coeffs(2) - z0 = surf%coeffs(3) - x1 = surf%coeffs(4) - y1 = surf%coeffs(5) - z1 = surf%coeffs(6) + x0 = surf % coeffs(1) + y0 = surf % coeffs(2) + z0 = surf % coeffs(3) + x1 = surf % coeffs(4) + y1 = surf % coeffs(5) + z1 = surf % coeffs(6) if (x >= x0 .and. x < x1 .and. y >= y0 .and. y < y1 .and. & & z >= z0 .and. z < z1) then s = SENSE_NEGATIVE @@ -672,8 +676,8 @@ contains c => cells(i) ! loop over each surface specification - do j = 1, c%n_items - index = c%boundary_list(j) + do j = 1, c % n_surfaces + index = c % surfaces(j) positive = (index > 0) index = abs(index) if (positive) then @@ -703,8 +707,8 @@ contains c => cells(i) ! loop over each surface specification - do j = 1, c%n_items - index = c%boundary_list(j) + do j = 1, c % n_surfaces + index = c % surfaces(j) positive = (index > 0) index = abs(index) diff --git a/src/global.f90 b/src/global.f90 index d9231ffe07..de41e590e1 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -54,8 +54,8 @@ module global type(ExtSource), target :: external_source ! Source and fission bank - type(Neutron), allocatable, target :: source_bank(:) - type(Bank), allocatable, target :: fission_bank(:) + type(Particle), allocatable, target :: source_bank(:) + type(Bank), allocatable, target :: fission_bank(:) ! Paths to input file, cross section data, etc character(100) :: & @@ -142,9 +142,9 @@ module global ! Particle type integer, parameter :: & - & NEUTRON_ = 1, & - & PHOTON_ = 2, & - & ELECTRON_ = 3 + & NEUTRON = 1, & + & PHOTON = 2, & + & ELECTRON = 3 ! Integer code for read error -- better hope this number is never ! used in an input file! diff --git a/src/main.f90 b/src/main.f90 index 19e966ff21..baedd5cd6a 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -90,7 +90,7 @@ contains integer :: i, j integer :: i_cycle integer :: i_particle - type(Neutron), pointer :: particle => null() + type(Particle), pointer :: p => null() CYCLE_LOOP: do i_cycle = 1, n_cycles @@ -99,8 +99,8 @@ contains HISTORY_LOOP: do j = 1, n_particles ! grab source particle from bank - particle => get_source_particle() - if ( .not. associated(particle) ) then + p => get_source_particle() + if ( .not. associated(p) ) then ! no particles left in source bank exit HISTORY_LOOP end if @@ -110,7 +110,7 @@ contains call RN_init_particle(int(i_particle,8)) ! transport particle - call transport(particle) + call transport(p) end do HISTORY_LOOP diff --git a/src/output.f90 b/src/output.f90 index dfc2f5b28b..4ec9f6b2a3 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -256,8 +256,8 @@ module output end if write(ou,*) ' Parent Cell = ' // int_to_str(c % parent) string = "" - do i = 1, c % n_items - select case (c % boundary_List(i)) + do i = 1, c % n_surfaces + select case (c % surfaces(i)) case (OP_LEFT_PAREN) string = trim(string) // ' (' case (OP_RIGHT_PAREN) @@ -267,7 +267,7 @@ module output case (OP_DIFFERENCE) string = trim(string) // ' !' case default - string = trim(string) // ' ' // int_to_str(c % boundary_list(i)) + string = trim(string) // ' ' // int_to_str(c % surfaces(i)) end select end do write(ou,*) ' Surface Specification:' // trim(string) diff --git a/src/physics.f90 b/src/physics.f90 index d9035de621..aed1a3cf07 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -2,7 +2,7 @@ module physics use global use geometry, only: find_cell, dist_to_boundary, cross_boundary - use types, only: Neutron + use types, only: Particle use mcnp_random, only: rang use output, only: error, message use search, only: binary_search @@ -17,9 +17,9 @@ contains ! geometry. !===================================================================== - subroutine transport(neut) + subroutine transport(p) - type(Neutron), pointer :: neut + type(Particle), pointer :: p character(250) :: msg integer :: i, surf @@ -31,46 +31,45 @@ contains real(8) :: distance real(8) :: Sigma ! total cross-section real(8) :: f ! interpolation factor - real(8) :: tmp(3) integer :: IE ! index on energy grid type(Universe), pointer :: univ - if (neut%cell == 0) then + if (p % cell == 0) then univ => universes(BASE_UNIVERSE) - call find_cell(univ, neut, found_cell) + call find_cell(univ, p, found_cell) - ! if neutron couldn't be located, print error + ! if particle couldn't be located, print error if (.not. found_cell) then - write(msg, 100) "Could not locate cell for neutron at: ", neut%xyz + write(msg, 100) "Could not locate cell for particle at: ", p % xyz 100 format (A,3ES11.3) call error(msg) end if end if if (verbosity >= 10) then - msg = "=== Particle " // trim(int_to_str(neut%uid)) // " ===" + msg = "=== Particle " // trim(int_to_str(p % uid)) // " ===" call message(msg, 10) - i = cells(neut%cell)%uid + 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 - neut%E = watt_spectrum(0.988_8, 2.249_8) + p % E = watt_spectrum(0.988_8, 2.249_8) ! find energy index, interpolation factor - call find_energy_index(neut) - IE = neut % IE - f = neut % interp + 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 (neut%alive) + do while (p % alive) ! Find the distance to the nearest boundary - call dist_to_boundary(neut, d_to_boundary, surf) + call dist_to_boundary(p, d_to_boundary, surf) ! Sample a distance to collision d_to_collision = -log(rang()) / 1.0 ! Sigma @@ -78,18 +77,18 @@ contains ! Select smaller of the two distances distance = min(d_to_boundary, d_to_collision) - ! Advance neutron - neut%xyz = neut%xyz + distance*neut%uvw + ! Advance particle + p%xyz = p%xyz + distance * p%uvw ! Add pathlength tallies if (d_to_collision > d_to_boundary) then - neut%surface = surf - neut%cell = 0 - call cross_boundary(neut) + p % surface = surf + p % cell = 0 + call cross_boundary(p) else ! collision - call collision(neut) + call collision(p) end if end do @@ -101,16 +100,16 @@ contains ! the interpolation factor for a particle at a certain energy !===================================================================== - subroutine find_energy_index(neut) + subroutine find_energy_index(p) - type(Neutron), pointer :: neut + type(Particle), pointer :: p integer :: IE real(8) :: E real(8) :: interp ! copy particle's energy - E = neut%E + E = p % E ! if particle's energy is outside of energy grid range, set to ! first or last index. Otherwise, do a binary search through the @@ -129,8 +128,8 @@ contains interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE)) ! set particle attributes - neut % IE = IE - neut % interp = interp + p % IE = IE + p % interp = interp end subroutine find_energy_index @@ -138,9 +137,9 @@ contains ! COLLISION !===================================================================== - subroutine collision(neut) + subroutine collision(p) - type(Neutron), pointer :: neut + type(Particle), pointer :: p type(AceContinuous), pointer :: table type(AceReaction), pointer :: rxn @@ -151,7 +150,7 @@ contains integer :: IE real(8) :: f, Sigma, total real(8) :: density, density_i - real(8) :: p + real(8) :: prob real(8), allocatable :: Sigma_t(:) ! tallies @@ -167,8 +166,8 @@ contains density_i = cMaterial%atom_percent(i)*density ! search nuclide energy grid - IE = table%grid_index(neut%IE) - f = (neut%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE)) + 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_t(i) = Sigma @@ -176,18 +175,18 @@ contains ! sample nuclide r1 = rang() - p = 0.0_8 + prob = 0.0_8 total = sum(Sigma_t) do i = 1, n_isotopes - p = p + Sigma_t(i) / total - if (r1 < p) exit + prob = prob + Sigma_t(i) / total + if (r1 < prob) exit end do ! Get table, total xs, interpolation factor table => xs_continuous(cMaterial%table(i)) Sigma = Sigma_t(i) - IE = table%grid_index(neut%IE) - f = (neut%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE)) + IE = table%grid_index(p % IE) + f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE)) density = cMaterial%atom_percent(i)*density ! free memory @@ -195,13 +194,14 @@ contains ! sample reaction channel r1 = rang()*Sigma - p = 0.0_8 + prob = 0.0_8 do i = 1, table%n_reaction rxn => table%reactions(i) if (rxn%MT >= 200) cycle if (IE < rxn%IE) cycle - p = p + density * (f*rxn%sigma(IE-rxn%IE+1) + (1-f)*(rxn%sigma(IE-rxn%IE+2))) - if (r1 < p) exit + prob = prob + density * (f*rxn%sigma(IE-rxn%IE+1) & + & + (1-f)*(rxn%sigma(IE-rxn%IE+2))) + if (r1 < prob) exit end do if (verbosity >= 10) then msg = " " // trim(reaction_name(rxn%MT)) // " with nuclide " // & @@ -212,11 +212,11 @@ contains ! call appropriate subroutine select case (rxn%MT) case (2) - call elastic_scatter(neut, table%awr) + call elastic_scatter(p, table%awr) case (102) - call n_gamma(neut) + call n_gamma(p) case default - call elastic_scatter(neut, table%awr) + call elastic_scatter(p, table%awr) end select end subroutine collision @@ -225,9 +225,9 @@ contains ! ELASTIC_SCATTER !===================================================================== - subroutine elastic_scatter(neut, awr) + subroutine elastic_scatter(p, awr) - type(Neutron), pointer :: neut + type(Particle), pointer :: p real(8), intent(in) :: awr real(8) :: phi ! azimuthal angle @@ -239,11 +239,11 @@ contains real(8) :: E integer :: IE - vel = sqrt(neut%E) + vel = sqrt(p % E) - vx = vel*neut%uvw(1) - vy = vel*neut%uvw(2) - vz = vel*neut%uvw(3) + vx = vel*p % uvw(1) + vy = vel*p % uvw(2) + vz = vel*p % uvw(3) vcx = vx/(awr + 1.0_8) @@ -277,13 +277,13 @@ contains E = vx*vx + vy*vy + vz*vz vel = sqrt(E) - neut%E = E - neut%uvw(1) = vx/vel - neut%uvw(2) = vy/vel - neut%uvw(3) = vz/vel + p % E = E + p % uvw(1) = vx/vel + p % uvw(2) = vy/vel + p % uvw(3) = vz/vel ! find energy index, interpolation factor - call find_energy_index(neut) + call find_energy_index(p) end subroutine elastic_scatter @@ -299,16 +299,16 @@ contains ! N_GAMMA !===================================================================== - subroutine n_gamma(neut) + subroutine n_gamma(p) - type(Neutron), pointer :: neut + type(Particle), pointer :: p integer :: cell_num character(250) :: msg - neut%alive = .false. + p % alive = .false. if (verbosity >= 10) then - cell_num = cells(neut%cell)%uid + cell_num = cells(p % cell)%uid msg = " Absorbed in cell " // trim(int_to_str(cell_num)) call message(msg, 10) end if diff --git a/src/source.f90 b/src/source.f90 index 1d73b45457..81b382f566 100644 --- a/src/source.f90 +++ b/src/source.f90 @@ -65,13 +65,13 @@ contains function get_source_particle() - type(Neutron), pointer :: get_source_particle + type(Particle), pointer :: get_source_particle ! increment index source_index = source_index + 1 ! if at end of bank, return null pointer - if ( source_index > size(source_bank) ) then + if (source_index > size(source_bank)) then get_source_particle => null() return end if diff --git a/src/types.f90 b/src/types.f90 index 6f83b5add6..14707d824e 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -7,14 +7,14 @@ module types !===================================================================== type Universe - integer :: uid - integer :: type - integer :: level - integer :: n_cells - integer, allocatable :: cells(:) - real(8) :: x0 - real(8) :: y0 - real(8) :: z0 + integer :: uid ! Unique ID + integer :: type ! Type + integer :: level ! Level of universe (0=base) + integer :: n_cells ! # of cells within + integer, allocatable :: cells(:) ! List of cells within + real(8) :: x0 ! Translation in x-coordinate + real(8) :: y0 ! Translation in y-coordinate + real(8) :: z0 ! Translation in z-coordinate end type Universe !===================================================================== @@ -22,10 +22,10 @@ module types !===================================================================== type Lattice - integer :: uid - integer :: type - integer :: level - real(8) :: pitch + integer :: uid ! Unique ID + integer :: type ! Type of lattice (square, hex, etc) + integer :: level ! Level of lattice + real(8) :: pitch ! Lattice pitch in cm end type Lattice !===================================================================== @@ -34,12 +34,13 @@ module types !===================================================================== type Surface - integer :: uid - integer :: type - real(8), allocatable :: coeffs(:) - integer, allocatable :: neighbor_pos(:) - integer, allocatable :: neighbor_neg(:) - integer :: bc + integer :: uid ! Unique ID + integer :: type ! Type of surface + real(8), allocatable :: coeffs(:) ! Definition of surface + integer, allocatable :: & + & neighbor_pos(:), & ! List of cells on positive side + & neighbor_neg(:) ! List of cells on negative side + integer :: bc ! Boundary condition end type Surface !===================================================================== @@ -47,23 +48,27 @@ module types !===================================================================== type Cell - integer :: uid - integer :: type - integer :: universe ! universe # this cell is in - integer :: fill ! universe # filling this cell - integer :: parent ! cell within which this cell resides - integer :: material - integer :: n_items - integer, allocatable :: boundary_list(:) + integer :: uid ! Unique ID + integer :: type ! Type of cell (normal, universe, lattice) + integer :: universe ! universe # this cell is in + integer :: fill ! universe # filling this cell + integer :: parent ! cell within which this cell resides + integer :: material ! Material within cell (0 for universe) + integer :: n_surfaces ! Number of surfaces within + integer, allocatable :: & + & surfaces(:) ! List of surfaces bounding cell -- note + ! that parentheses, union, etc operators + ! will be listed here too end type Cell !===================================================================== -! NEUTRON describes the state of a neutron being transported through +! PARTICLE describes the state of a particle being transported through ! the geometry !===================================================================== - type Neutron + type Particle integer :: uid ! Unique ID + integer :: type ! Particle type (n, p, e, etc) real(8) :: xyz(3) ! location real(8) :: uvw(3) ! directional cosines real(8) :: E ! energy @@ -74,7 +79,7 @@ module types integer :: surface ! current surface real(8) :: wgt ! particle weight logical :: alive ! is particle alive? - end type Neutron + end type Particle !===================================================================== ! BANK is used for storing fission sites in criticality