diff --git a/ChangeLog b/ChangeLog index 7f6e23350..249a19c34 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,19 @@ +2011-03-28 Paul Romano + + * fileio.f90: Ability to read lattice entries. + * geometry.f90: Added lattice capabilities. Changed cross_boundary + to cross_surface. Added many comments. In many places, changed + testing of c%fill to c%type instead since fill could be universe + or lattice. + * Makefile: Added -traceback + * output.f90: Added print_lattice, other minor changes. (Fixed + indentation so looks like everything changed) + * physics.f90: Call cross_lattice if crossing lattice boundary + instead of surface. + * source.f90: Set values for xyz_local, lattice indices + * types.f90: Extending Lattice type, modified Particle type to + include local coordinates, current lattice, and lattice indices. + 2011-03-25 Paul Romano * physics.f90: Resample outgoing energy if sampled energy is diff --git a/src/Makefile b/src/Makefile index b9e0b2e16..0737d678d 100644 --- a/src/Makefile +++ b/src/Makefile @@ -26,7 +26,7 @@ test_objects = $(modules:.f90=.o) $(test:.f90=.o) #-------------------------------------------------------------------- # Compiler Options F90 = ifort -F90FLAGS = -g -fpp +F90FLAGS = -g -fpp -traceback MPI = /opt/mpich2/1.3.3-gcc MPIF90 = $(MPI)/bin/mpif90 MPI_COMPILE_FLAGS = $(shell mpif90 -f90='' -compile_info) diff --git a/src/fileio.f90 b/src/fileio.f90 index ba4b20b1e..65ef02318 100644 --- a/src/fileio.f90 +++ b/src/fileio.f90 @@ -97,6 +97,7 @@ contains call dict_create(ucount_dict) call dict_create(universe_dict) + call dict_create(lattice_dict) open(FILE=filename, UNIT=in, STATUS='old', & & ACTION='read', IOSTAT=ioError) @@ -131,6 +132,15 @@ contains n_surfaces = n_surfaces + 1 case ('material') n_materials = n_materials + 1 + case ('lattice') + n_lattices = n_lattices + 1 + + ! For lattices, 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(2)) + call dict_add_key(lattice_dict, universe_num, n_lattices) + end select end do @@ -155,6 +165,7 @@ contains ! Allocate arrays for cells, surfaces, etc. allocate(cells(n_cells)) allocate(universes(n_universes)) + allocate(lattices(n_lattices)) allocate(surfaces(n_surfaces)) allocate(materials(n_materials)) @@ -251,7 +262,7 @@ contains call read_cell(index_cell, words, n) case ('surface') - ! Read data + ! Read data from surface entry index_surface = index_surface + 1 call read_surface(index_surface, words, n) @@ -261,6 +272,7 @@ contains call read_lattice(index_lattice, words, n) case ('material') + ! Read data from material entry index_material = index_material + 1 call read_material(index_material, words, n) @@ -329,7 +341,7 @@ contains end do ! adjust material indices - if (c % fill == 0) then + if (c % type == CELL_NORMAL) then index = dict_get_key(material_dict, c%material) if (index == DICT_NULL) then msg = "Could not find material " // trim(int_to_str(c%material)) // & @@ -354,9 +366,13 @@ contains integer, intent(in) :: level ! level of universe integer :: i ! index for cells in universe + integer :: x,y ! indices for lattice positions integer :: index ! index in cells array + integer :: universe_num type(Cell), pointer :: c => null() type(Universe), pointer :: subuniverse => null() + type(Lattice), pointer :: lat => null() + type(DictionaryII), pointer :: dict => null() ! set level of the universe univ % level = level @@ -369,10 +385,31 @@ contains ! if this cell is filled with another universe, recursively ! call this subroutine - if (c % fill > 0) then + if (c % type == CELL_FILL) then subuniverse => universes(c % fill) call build_universe(subuniverse, index, level + 1) end if + + ! if this cell is filled by a lattice, need to build the + ! universe for each unique lattice element + if (c % type == CELL_LATTICE) then + lat => lattices(c % fill) + call dict_create(dict) + do x = 1, lat % n_x + do y = 1, lat % n_y + universe_num = lat % element(x,y) + if (dict_has_key(dict, universe_num)) then + cycle + else + call dict_add_key(dict, universe_num, 0) + + subuniverse => universes(universe_num) + call build_universe(subuniverse, index, level + 1) + end if + end do + end do + end if + end do end subroutine build_universe @@ -415,7 +452,6 @@ contains ! Read cell material if (trim(words(4)) == 'fill') then - c % type = CELL_FILL c % material = 0 n_surfaces = n_words - 5 @@ -425,8 +461,15 @@ contains msg = "Invalid universe fill: " // words(5) call error(msg) end if - c % fill = dict_get_key(universe_dict, universe_num) + ! check whether universe or lattice + if (dict_has_key(universe_dict, universe_num)) then + c % type = CELL_FILL + c % fill = dict_get_key(universe_dict, universe_num) + elseif (dict_has_key(lattice_dict, universe_num)) then + c % type = CELL_LATTICE + c % fill = dict_get_key(lattice_dict, universe_num) + end if else c % type = CELL_NORMAL c % material = str_to_int(words(4)) @@ -439,7 +482,7 @@ contains end if ! Assign number of items - c%n_surfaces = n_surfaces + c % n_surfaces = n_surfaces ! Read list of surfaces allocate(c%surfaces(n_surfaces)) @@ -560,9 +603,83 @@ contains subroutine read_lattice(index, words, n_words) - integer, intent(in) :: index - character(*), intent(in) :: words(n_words) - integer, intent(in) :: n_words + integer, intent(in) :: index ! index in lattices array + character(*), intent(in) :: words(n_words) ! words in lattice entry + integer, intent(in) :: n_words ! number of words + + integer :: universe_num ! user-specified universe number + integer :: n_x ! number of lattice cells in x direction + integer :: n_y ! number of lattice cells in y direction + integer :: i,j ! loop indices for Lattice % universes + integer :: index_word ! index in words array + character(250) :: msg ! output/error/message + character(32) :: word ! single word + type(Lattice), pointer :: lat => null() + + lat => lattices(index) + + ! Read lattice universe + universe_num = str_to_int(words(2)) + if (universe_num == ERROR_CODE) then + msg = "Invalid universe: " // words(2) + call error(msg) + end if + lat % uid = universe_num + + ! Read lattice type + word = words(3) + call lower_case(word) + select case(trim(word)) + case ('rect') + lat % type = LATTICE_RECT + case ('hex') + lat % type = LATTICE_HEX + case default + msg = "Invalid lattice type: " // words(3) + call error(msg) + end select + + ! 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 + msg = "Invalid number of lattice cells in x-direction: " // words(4) + call error(msg) + elseif (n_y == ERROR_CODE) then + msg = "Invalid number of lattice cells in y-direction: " // words(5) + call error(msg) + end if + lat % n_x = n_x + lat % n_y = n_y + + ! Read lattice origin coordinates + lat % x0 = str_to_real(words(6)) + lat % y0 = str_to_real(words(7)) + + ! Read lattice cell widths + lat % width_x = str_to_real(words(8)) + lat % width_y = str_to_real(words(9)) + + ! Make sure enough lattice positions specified + if (n_words - 9 < n_x * n_y) then + print *, n_words, n_x, n_y + msg = "Not enough lattice positions specified." + call error(msg) + end if + + ! Read lattice positions + allocate(lat % element(n_x,n_y)) + do j = 0, n_y - 1 + 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 + msg = "Invalid universe number: " // words(index_word) + call error(msg) + end if + lat % element(i, n_y-j) = dict_get_key(universe_dict, universe_num) + end do + end do end subroutine read_lattice diff --git a/src/geometry.f90 b/src/geometry.f90 index 761464e53..b7140a820 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -19,15 +19,15 @@ contains type(Particle), pointer :: p logical :: in_cell - integer, allocatable :: expression(:) - integer :: specified_sense - integer :: actual_sense - integer :: n_surfaces - integer :: i, j - integer :: surf_num - integer :: current_surface + integer, allocatable :: expression(:) ! copy of surfaces list + integer :: specified_sense ! specified sense of surface in list + integer :: actual_sense ! sense of particle wrt surface + integer :: n_surfaces ! number of surfaces in cell + integer :: i ! index of surfaces in cell + integer :: surf_num ! index in surfaces array (with sign) + integer :: current_surface ! current surface of particle (with sign) + character(250) :: msg ! output/error message type(Surface), pointer :: surf => null() - character(250) :: msg current_surface = p%surface @@ -59,7 +59,7 @@ contains ! Compare sense of point to specified sense specified_sense = sign(1,expression(i)) - call sense(surf, p%xyz, actual_sense) + actual_sense = sense(surf, p%xyz_local) if (actual_sense == specified_sense) then expression(i) = 1 else @@ -95,7 +95,9 @@ contains character(250) :: msg ! error message integer :: i ! index over cells + integer :: x, y type(Cell), pointer :: c ! pointer to cell + type(Lattice), pointer :: lat ! pointer to lattice type(Universe), pointer :: lower_univ ! if particle is in lower ! universe, use this pointer ! to call recursively @@ -107,9 +109,20 @@ contains c => cells(univ % cells(i)) if (cell_contains(c, p)) then - ! If this cell contains a universe of lattice, search for + ! If this cell contains a universe or lattice, search for ! the particle in that universe/lattice - if (c % fill > 0) then + if (c % type == CELL_NORMAL) then + ! set current pointers + found = .true. + cCell => c + cMaterial => materials(cCell%material) + cUniverse => univ + + ! set particle attributes + p%cell = univ % cells(i) + p%universe = dict_get_key(universe_dict, univ % uid) + exit + elseif (c % type == CELL_FILL) then lower_univ => universes(c % fill) call find_cell(lower_univ, p, found) if (found) then @@ -118,17 +131,34 @@ contains msg = "Could not locate particle in universe: " call error(msg) end if - else - ! set current pointers - found = .true. - cCell => c - cMaterial => materials(cCell%material) - cUniverse => univ + elseif (c % type == CELL_LATTICE) then + ! Set current lattice + lat => lattices(c % fill) + cLattice => lat + p % lattice = c % fill + + ! determine universe based on lattice position + x = ceiling((p%xyz(1) - lat%x0)/lat%width_x) + y = ceiling((p%xyz(2) - lat%y0)/lat%width_y) + lower_univ => universes(lat % element(x,y)) + + ! adjust local position of particle + p%xyz_local(1) = p%xyz(1) - (lat%x0 + (x-0.5)*lat%width_x) + p%xyz_local(2) = p%xyz(2) - (lat%y0 + (y-0.5)*lat%width_y) + p%xyz_local(3) = p%xyz(3) - ! set particle attributes - p%cell = dict_get_key(cell_dict, c % uid) - p%universe = dict_get_key(universe_dict, univ % uid) - exit + ! set particle lattice indices + p % index_x = x + p % index_y = y + + call find_cell(lower_univ, p, found) + if (found) then + exit + else + msg = "Could not locate particle in lattice: " & + & // int_to_str(lat % uid) + call error(msg) + end if end if end if @@ -137,19 +167,19 @@ contains end subroutine find_cell !===================================================================== -! CROSS_BOUNDARY moves a particle into a new cell +! CROSS_SURFACE moves a particle into a new cell !===================================================================== - subroutine cross_boundary(p) + subroutine cross_surface(p) type(Particle), pointer :: p - type(Surface), pointer :: surf - type(Cell), pointer :: c - integer :: i - integer :: index_cell - character(250) :: msg - logical :: found + integer :: i ! index of neighbors + integer :: index_cell ! index in cells array + logical :: found ! particle found in universe? + character(250) :: msg ! output/error message? + type(Surface), pointer :: surf + type(Cell), pointer :: c type(Universe), pointer :: lower_univ => null() surf => surfaces(abs(p%surface)) @@ -175,7 +205,7 @@ contains index_cell = surf%neighbor_pos(i) c => cells(index_cell) if (cell_contains(c, p)) then - if (c % fill > 0) then + if (c % type == CELL_FILL) then lower_univ => universes(c % fill) call find_cell(lower_univ, p, found) if (.not. found) then @@ -198,7 +228,7 @@ contains index_cell = surf%neighbor_neg(i) c => cells(index_cell) if (cell_contains(c, p)) then - if (c % fill > 0) then + if (c % type == CELL_FILL) then lower_univ => universes(c % fill) call find_cell(lower_univ, p, found) if (.not. found) then @@ -233,44 +263,158 @@ contains & ", it could not be located in any cell and it did not leak." call error(msg) - end subroutine cross_boundary + end subroutine cross_surface !===================================================================== -! DIST_TO_BOUNDARY calculates the distance to the nearest boundary of -! 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. +! CROSS_LATTICE moves a particle into a new lattice element !===================================================================== - subroutine dist_to_boundary(p, dist, surf, other_cell) + subroutine cross_lattice(p) - type(Particle), pointer :: p - real(8), intent(out) :: dist - integer, intent(out) :: surf - integer, optional, intent(in) :: other_cell + type(Particle), pointer :: p + integer :: i_x ! x index in lattice + integer :: i_y ! y index in lattice + real(8) :: d_left ! distance to left side + real(8) :: d_right ! distance to right side + real(8) :: d_bottom ! distance to bottom side + real(8) :: d_top ! distance to top side + real(8) :: dist ! shortest distance + real(8) :: x ! x coordinate in local lattice element + real(8) :: y ! y coordinate in local lattice element + real(8) :: z ! z coordinate in local lattice element + real(8) :: u ! cosine of angle with x axis + real(8) :: v ! cosine of angle with y axis + real(8) :: x0 ! half the width of lattice element + real(8) :: y0 ! half the height of lattice element + logical :: found ! particle found in cell? + character(250) :: msg ! output/error message + type(Lattice), pointer :: lat + type(Universe), pointer :: univ + + if (verbosity >= 10) then + msg = " Crossing lattice" + call message(msg, 10) + end if + + lat => lattices(p % lattice) + + u = p % uvw(1) + v = p % uvw(2) + + if (lat % type == LATTICE_RECT) then + x = p % xyz_local(1) + y = p % xyz_local(2) + z = p % xyz_local(3) + x0 = lat % width_x * 0.5 + y0 = lat % width_y * 0.5 + + dist = INFINITY + + ! left and right sides + if (u > 0) then + d_left = INFINITY + d_right = (x0 - x)/u + else + d_left = -(x0 + x)/u + d_right = INFINITY + end if + + ! top and bottom sides + if (v > 0) then + d_bottom = INFINITY + d_top = (y0 - y)/v + else + d_bottom = -(y0 + y)/v + d_top = INFINITY + end if + + dist = min(d_left, d_right, d_top, d_bottom) + if (dist == d_left) then + ! Move particle to left element + p % index_x = p % index_x - 1 + p % xyz_local(1) = x0 + + elseif (dist == d_right) then + ! Move particle to right element + p % index_x = p % index_x + 1 + p % xyz_local(1) = -x0 + + elseif (dist == d_bottom) then + ! Move particle to bottom element + p % index_y = p % index_y - 1 + p % xyz_local(2) = y0 + + elseif (dist == d_top) then + ! Move particle to top element + p % index_y = p % index_y + 1 + p % xyz_local(2) = -y0 + + end if + elseif (lat % type == LATTICE_HEX) then + ! TODO: Add hex lattice support + end if + + ! Check to make sure still in lattice + i_x = p % index_x + i_y = p % index_y + if (i_x < 1 .or. i_x > lat % n_x) then + msg = "Reached edge of lattice." + call error(msg) + elseif (i_y < 1 .or. i_y > lat % n_y) then + msg = "Reached edge of lattice." + call error(msg) + end if + + ! Find universe for next lattice element + univ => universes(lat % element(i_x,i_y)) + + ! Find cell in next lattice element + call find_cell(univ, p, found) + if (.not. found) then + msg = "Could not locate particle in universe: " + call error(msg) + end if + + end subroutine cross_lattice + +!===================================================================== +! DIST_TO_BOUNDARY calculates the distance to the nearest boundary 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(p, dist, surf, in_lattice) + + type(Particle), pointer :: p + real(8), intent(out) :: dist + integer, intent(out) :: surf + logical, intent(out) :: in_lattice + + integer, allocatable :: expression(:) ! copy of surface list + integer :: i ! index for surface in cell + integer :: n_surfaces ! total number of surfaces to check + integer :: n1 ! number of surfaces in current cell + integer :: n2 ! number of surfaces in parent cell + integer :: index_surf ! index in surfaces array (with sign) + integer :: current_surf ! current surface + real(8) :: x,y,z ! particle coordinates + real(8) :: u,v,w ! particle directions + real(8) :: d ! evaluated distance + real(8) :: x0,y0,z0 ! coefficients for surface + real(8) :: r ! radius for quadratic surfaces + real(8) :: tmp ! dot product of surface normal with direction + real(8) :: a,b,c,k ! quadratic equation coefficients + real(8) :: quad ! discriminant of quadratic equation + logical :: on_surface ! is particle on surface? + character(250) :: msg ! output/error message type(Cell), pointer :: cell_p => null() type(Cell), pointer :: parent_p => null() type(Surface), pointer :: surf_p => null() - integer :: i - integer :: n_surfaces, n1, n2 - integer, allocatable :: expression(:) - integer :: index_surf - integer :: current_surf - real(8) :: x,y,z,u,v,w - real(8) :: d - real(8) :: x0,y0,z0,r - real(8) :: tmp - real(8) :: a,b,c,k - real(8) :: quad - character(250) :: msg - logical :: on_surface + type(LatticE), pointer :: lat => null() - if (present(other_cell)) then - cell_p => cells(other_cell) - else - cell_p => cells(p%cell) - end if + cell_p => cells(p%cell) current_surf = p%surface @@ -290,16 +434,24 @@ contains expression(n1+1:n1+n2) = parent_p % surfaces end if + u = p % uvw(1) + v = p % uvw(2) + w = p % uvw(3) + ! loop over all surfaces dist = INFINITY 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) - + if (i <= n1) then + ! in local cell, so use xyz_local + x = p % xyz_local(1) + y = p % xyz_local(2) + z = p % xyz_local(3) + else + ! in parent cell, so use xyz + x = p % xyz(1) + y = p % xyz(2) + z = p % xyz(3) + end if ! check for coincident surface index_surf = expression(i) if (index_surf == current_surf) then @@ -570,6 +722,45 @@ contains end do + ! Check lattice surfaces + in_lattice = .false. + if (p % lattice > 0) then + lat => lattices(p % lattice) + if (lat % type == LATTICE_RECT) then + x = p % xyz_local(1) + y = p % xyz_local(2) + z = p % xyz_local(3) + x0 = lat % width_x * 0.5 + y0 = lat % width_y * 0.5 + + + ! left and right sides + if (u > 0) then + d = (x0 - x)/u + else + d = -(x + x0)/u + end if + if (d < dist) then + dist = d + in_lattice = .true. + end if + + ! top and bottom sides + if (v > 0) then + d = (y0 - y)/v + else + d = -(y + y0)/v + end if + if (d < dist) then + dist = d + in_lattice = .true. + end if + + elseif (lat % type == LATTICE_HEX) then + ! TODO: Add hex lattice support + end if + end if + ! deallocate expression deallocate(expression) @@ -581,17 +772,27 @@ contains ! a particular point is in. !===================================================================== - subroutine sense(surf, xyz, s) + function sense(surf, xyz) result(s) - type(Surface), intent(in) :: surf - real(8), intent(in) :: xyz(3) - integer, intent(out) :: s + type(Surface), pointer :: surf ! surface + real(8), intent(in) :: xyz(3) ! coordinates of particle + integer :: s ! sense of particle - real(8) :: x,y,z - real(8) :: func - real(8) :: A,B,C,D,E,F,G,H,I,J - real(8) :: x0, y0, z0, r - real(8) :: x1, y1, z1 + real(8) :: x,y,z ! coordinates of particle + real(8) :: func ! surface function evaluated at point + real(8) :: A ! coefficient on x**2 term in GQ + real(8) :: B ! coefficient on y**2 term in GQ + real(8) :: C ! coefficient on z**2 term in GQ + real(8) :: D ! coefficient on x*y term in GQ + real(8) :: E ! coefficient on y*z term in GQ + real(8) :: F ! coefficient on x*z term in GQ + real(8) :: G ! coefficient on x term in GQ + real(8) :: H ! coefficient on y term in GQ + real(8) :: I ! coefficient on z term in GQ + real(8) :: J ! coefficient on constant term in GQ + real(8) :: x0,y0,z0 ! coefficients for quadratic surfaces / box + real(8) :: r ! radius for quadratic surfaces + real(8) :: x1,y1,z1 ! upper-right corner of box x = xyz(1) y = xyz(2) @@ -715,7 +916,7 @@ contains s = SENSE_NEGATIVE end if - end subroutine sense + end function sense !===================================================================== ! NEIGHBOR_LISTS builds a list of neighboring cells to each surface to @@ -724,16 +925,15 @@ contains subroutine neighbor_lists() + integer :: i ! index in cells/surfaces array + integer :: j ! index of surface in cell + integer :: index ! index in count arrays + integer, allocatable :: count_positive(:) ! # of cells on positive side + integer, allocatable :: count_negative(:) ! # of cells on negative side + logical :: positive ! positive side specified in surface list + character(250) :: msg ! output/error message type(Cell), pointer :: c type(Surface), pointer :: surf - integer :: i, j - integer :: index - integer :: surf_num - logical :: positive - character(250) :: msg - - integer, allocatable :: count_positive(:) - integer, allocatable :: count_negative(:) msg = "Building neighboring cells lists for each surface..." call message(msg, 4) diff --git a/src/main.f90 b/src/main.f90 index 336433c23..5cda3dc7f 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -5,7 +5,7 @@ program main & normalize_ao, build_universe use output, only: title, echo_input, message, warning, error, & & print_summary, print_particle - use geometry, only: sense, cell_contains, neighbor_lists + use geometry, only: 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 @@ -81,7 +81,10 @@ program main call init_source() ! start problem +! surfaces(1)%bc = BC_VACUUM +! surfaces(2)%bc = BC_VACUUM surfaces(1)%bc = BC_VACUUM +! surfaces(4)%bc = BC_VACUUM call run_problem() ! deallocate arrays diff --git a/src/output.f90 b/src/output.f90 index 0aa2a4ee2..7d4c2addf 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -11,7 +11,7 @@ module output integer :: ou = OUTPUT_UNIT integer :: eu = ERROR_UNIT - contains +contains !===================================================================== ! TITLE prints the main title banner as well as information about the @@ -19,131 +19,129 @@ module output ! run. !===================================================================== - subroutine title() + subroutine title() - character(10) :: date - character(8) :: time + character(10) :: date + character(8) :: time - write(ou,*) - write(ou,*) ' .d88888b. 888b d888 .d8888b. ' - write(ou,*) 'd88P" "Y88b 8888b d8888 d88P Y88b' - write(ou,*) '888 888 88888b.d88888 888 888' - write(ou,*) '888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ' - write(ou,*) '888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ' - write(ou,*) '888 888 888 888 88888888 888 888 888 Y8P 888 888 888' - write(ou,*) 'Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P' - write(ou,*) ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P" ' - write(ou,*) '____________888________________________________________________' - write(ou,*) ' 888 ' - write(ou,*) ' 888 ' - write(ou,*) + write(ou,*) + write(ou,*) ' .d88888b. 888b d888 .d8888b. ' + write(ou,*) 'd88P" "Y88b 8888b d8888 d88P Y88b' + write(ou,*) '888 888 88888b.d88888 888 888' + write(ou,*) '888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ' + write(ou,*) '888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ' + write(ou,*) '888 888 888 888 88888888 888 888 888 Y8P 888 888 888' + write(ou,*) 'Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P' + write(ou,*) ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P" ' + write(ou,*) '____________888________________________________________________' + write(ou,*) ' 888 ' + write(ou,*) ' 888 ' + write(ou,*) - ! Write version information -! write(ou,*) ' ______________________________________________________' - write(ou,*) ' Developed At: Massachusetts Institute of Technology' - write(ou,*) ' Lead Developer: Paul K. Romano' - write(ou,100) VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE -100 format (6X,"Version:",9X,I1,".",I1,".",I1) + ! Write version information + write(ou,*) ' Developed At: Massachusetts Institute of Technology' + write(ou,*) ' Lead Developer: Paul K. Romano' + write(ou,100) VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE +100 format (6X,"Version:",9X,I1,".",I1,".",I1) - ! Write the date and time - call get_today( date, time ) - write(ou,101) trim(date), trim(time) -101 format (6X,"Date/Time:",7X,A,1X,A) -! write(ou,*) ' ______________________________________________________' - write(ou,*) + ! Write the date and time + call get_today( date, time ) + write(ou,101) trim(date), trim(time) +101 format (6X,"Date/Time:",7X,A,1X,A) + write(ou,*) - end subroutine title + end subroutine title !===================================================================== ! ECHO_INPUT displays summary information about the problem about to ! be run after reading all the input. !===================================================================== - subroutine echo_input() + subroutine echo_input() - character(32) :: string + character(32) :: string - ! Display problem summary - write(ou,*) '=============================================' - write(ou,*) '=> PROBLEM SUMMARY <=' - write(ou,*) '=============================================' - write(ou,*) - if ( problem_type == PROB_CRITICALITY ) then - write(ou,100) 'Problem type:', 'Criticality' - write(ou,100) 'Number of Cycles:', int_to_str(n_cycles) - write(ou,100) 'Number of Inactive Cycles:', int_to_str(n_inactive) - elseif ( problem_type == PROB_SOURCE ) then - write(ou,100) 'Problem type:', 'External Source' - end if - write(ou,100) 'Number of Particles:', int_to_str(n_particles) - write(ou,*) - ! Display geometry summary - write(ou,*) '=============================================' - write(ou,*) '=> GEOMETRY SUMMARY <=' - write(ou,*) '=============================================' - write(ou,*) - write(ou,100) 'Number of Cells:', int_to_str(n_cells) - write(ou,100) 'Number of Surfaces:', int_to_str(n_surfaces) - write(ou,100) 'Number of Materials:', int_to_str(n_materials) - write(ou,*) + ! Display problem summary + write(ou,*) '=============================================' + write(ou,*) '=> PROBLEM SUMMARY <=' + write(ou,*) '=============================================' + write(ou,*) + if ( problem_type == PROB_CRITICALITY ) then + write(ou,100) 'Problem type:', 'Criticality' + write(ou,100) 'Number of Cycles:', int_to_str(n_cycles) + write(ou,100) 'Number of Inactive Cycles:', int_to_str(n_inactive) + elseif ( problem_type == PROB_SOURCE ) then + write(ou,100) 'Problem type:', 'External Source' + end if + write(ou,100) 'Number of Particles:', int_to_str(n_particles) + write(ou,*) + ! Display geometry summary + write(ou,*) '=============================================' + write(ou,*) '=> GEOMETRY SUMMARY <=' + write(ou,*) '=============================================' + write(ou,*) + write(ou,100) 'Number of Cells:', int_to_str(n_cells) + write(ou,100) 'Number of Surfaces:', int_to_str(n_surfaces) + write(ou,100) 'Number of Materials:', int_to_str(n_materials) + write(ou,*) - ! Format descriptor for columns -100 format (1X,A,T35,A) + ! Format descriptor for columns +100 format (1X,A,T35,A) - end subroutine echo_input + end subroutine echo_input !===================================================================== ! MESSAGE displays an informational message to the log file and the ! standard output stream. !===================================================================== - subroutine message( msg, level ) + subroutine message( msg, level ) - character(*), intent(in) :: msg - integer, intent(in) :: level + character(*), intent(in) :: msg + integer, intent(in) :: level - integer :: n_lines - integer :: i + integer :: n_lines + integer :: i - ! Only allow master to print to screen - if (.not. master) return + ! Only allow master to print to screen + if (.not. master) return - if ( level <= verbosity ) then - n_lines = (len_trim(msg)-1)/79 + 1 - do i = 1, n_lines - write(ou, fmt='(1X,A79)') msg(79*(i-1)+1:79*i) - end do - end if + if ( level <= verbosity ) then + n_lines = (len_trim(msg)-1)/79 + 1 + do i = 1, n_lines + write(ou, fmt='(1X,A79)') msg(79*(i-1)+1:79*i) + end do + end if - end subroutine message + end subroutine message !===================================================================== ! WARNING issues a warning to the user in the log file and the ! standard output stream. !===================================================================== - subroutine warning( msg ) + subroutine warning( msg ) - character(*), intent(in) :: msg + character(*), intent(in) :: msg - integer :: n_lines - integer :: i + integer :: n_lines + integer :: i - ! Only allow master to print to screen - if (.not. master) return + ! Only allow master to print to screen + if (.not. master) return - write(ou, fmt='(1X,A9)', advance='no') 'WARNING: ' + write(ou, fmt='(1X,A9)', advance='no') 'WARNING: ' - n_lines = (len_trim(msg)-1)/70 + 1 - do i = 1, n_lines - if ( i == 1 ) then - write(ou, fmt='(A70)') msg(70*(i-1)+1:70*i) - else - write(ou, fmt='(10X,A70)') msg(70*(i-1)+1:70*i) - end if - end do - - end subroutine warning + n_lines = (len_trim(msg)-1)/70 + 1 + do i = 1, n_lines + if ( i == 1 ) then + write(ou, fmt='(A70)') msg(70*(i-1)+1:70*i) + else + write(ou, fmt='(10X,A70)') msg(70*(i-1)+1:70*i) + end if + end do + + end subroutine warning !===================================================================== ! ERROR alerts the user that an error has been encountered and @@ -151,336 +149,369 @@ module output ! considered 'fatal' and hence the program is aborted. !===================================================================== - subroutine error( msg ) + subroutine error( msg ) - character(*), intent(in) :: msg + character(*), intent(in) :: msg - integer :: n_lines - integer :: i + integer :: n_lines + integer :: i - ! Only allow master to print to screen - if (master) then - write(eu, fmt='(1X,A7)', advance='no') 'ERROR: ' + ! Only allow master to print to screen + if (master) then + write(eu, fmt='(1X,A7)', advance='no') 'ERROR: ' - n_lines = (len_trim(msg)-1)/72 + 1 - do i = 1, n_lines - if ( i == 1 ) then - write(eu, fmt='(A72)') msg(72*(i-1)+1:72*i) - else - write(eu, fmt='(7X,A72)') msg(72*(i-1)+1:72*i) - end if - end do - write(eu,*) - end if + n_lines = (len_trim(msg)-1)/72 + 1 + do i = 1, n_lines + if ( i == 1 ) then + write(eu, fmt='(A72)') msg(72*(i-1)+1:72*i) + else + write(eu, fmt='(7X,A72)') msg(72*(i-1)+1:72*i) + end if + end do + write(eu,*) + end if - ! All processors abort - call free_memory() - - end subroutine error + ! All processors abort + call free_memory() + + end subroutine error !===================================================================== ! GET_TODAY determines the date and time at which the program began ! execution and returns it in a readable format !===================================================================== - subroutine get_today(today_date, today_time) + subroutine get_today(today_date, today_time) - character(10), intent(out) :: today_date - character(8), intent(out) :: today_time + character(10), intent(out) :: today_date + character(8), intent(out) :: today_time - character(8) :: date - character(10) :: time - character(5) :: zone - integer :: val(8) + character(8) :: date + character(10) :: time + character(5) :: zone + integer :: val(8) - call date_and_time(date, time, zone, val) - ! val(1) = year (YYYY) - ! val(2) = month (MM) - ! val(3) = day (DD) - ! val(4) = timezone - ! val(5) = hours (HH) - ! val(6) = minutes (MM) - ! val(7) = seconds (SS) - ! val(8) = milliseconds - - if ( val(2) < 10 ) then - if ( val(3) < 10 ) then - today_date = date(6:6) // "/" // date(8:8) // "/" // date(1:4) - else - today_date = date(6:6) // "/" // date(7:8) // "/" // date(1:4) - end if - else - if ( val(3) < 10 ) then - today_date = date(5:6) // "/" // date(8:8) // "/" // date(1:4) - else - today_date = date(5:6) // "/" // date(7:8) // "/" // date(1:4) - end if - end if - today_time = time(1:2) // ":" // time(3:4) // ":" // time(5:6) + call date_and_time(date, time, zone, val) + ! val(1) = year (YYYY) + ! val(2) = month (MM) + ! val(3) = day (DD) + ! val(4) = timezone + ! val(5) = hours (HH) + ! val(6) = minutes (MM) + ! val(7) = seconds (SS) + ! val(8) = milliseconds - end subroutine get_today + if ( val(2) < 10 ) then + if ( val(3) < 10 ) then + today_date = date(6:6) // "/" // date(8:8) // "/" // date(1:4) + else + today_date = date(6:6) // "/" // date(7:8) // "/" // date(1:4) + end if + else + if ( val(3) < 10 ) then + today_date = date(5:6) // "/" // date(8:8) // "/" // date(1:4) + else + today_date = date(5:6) // "/" // date(7:8) // "/" // date(1:4) + end if + end if + today_time = time(1:2) // ":" // time(3:4) // ":" // time(5:6) + + end subroutine get_today !===================================================================== ! PRINT_PARTICLE displays the attributes of a particle !===================================================================== - subroutine print_particle(p) + subroutine print_particle(p) - type(Particle), pointer :: p + type(Particle), pointer :: p - integer :: i - type(Cell), pointer :: c => null() - type(Surface), pointer :: s => null() - type(Universe), pointer :: u => null() - character(250) :: string + integer :: i + type(Cell), pointer :: c => null() + type(Surface), pointer :: s => null() + type(Universe), pointer :: u => null() + character(250) :: string - select case (p % type) - case (NEUTRON) - write(ou,*) 'Neutron ' // int_to_str(p % uid) - case (PHOTON) - write(ou,*) 'Photon ' // int_to_str(p % uid) - case (ELECTRON) - write(ou,*) 'Electron ' // int_to_str(p % uid) - case default - write(ou,*) 'Unknown Particle ' // int_to_str(p % uid) - end select - write(ou,100) 'x = ', p % xyz(1) - write(ou,100) 'y = ', p % xyz(2) - write(ou,100) 'z = ', p % xyz(3) - write(ou,100) 'u = ', p % uvw(1) - write(ou,100) 'v = ', p % uvw(2) - write(ou,100) 'w = ', p % uvw(3) - write(ou,100) 'Weight = ', p % wgt - write(ou,100) 'Energy = ', p % E - write(ou,*) ' IE = ' // int_to_str(p % IE) - write(ou,100) 'Interpolation factor = ', p % interp - - if (p % cell > 0) then - c => cells(p % cell) - write(ou,*) ' Cell = ' // int_to_str(c % uid) - else - write(ou,*) ' Cell not determined' - end if + select case (p % type) + case (NEUTRON) + write(ou,*) 'Neutron ' // int_to_str(p % uid) + case (PHOTON) + write(ou,*) 'Photon ' // int_to_str(p % uid) + case (ELECTRON) + write(ou,*) 'Electron ' // int_to_str(p % uid) + case default + write(ou,*) 'Unknown Particle ' // int_to_str(p % uid) + end select + write(ou,100) 'x = ', p % xyz(1) + write(ou,100) 'y = ', p % xyz(2) + write(ou,100) 'z = ', p % xyz(3) + write(ou,100) 'x local = ', p % xyz_local(1) + write(ou,100) 'y local = ', p % xyz_local(2) + write(ou,100) 'z local = ', p % xyz_local(3) + write(ou,100) 'u = ', p % uvw(1) + write(ou,100) 'v = ', p % uvw(2) + write(ou,100) 'w = ', p % uvw(3) + write(ou,100) 'Weight = ', p % wgt + write(ou,100) 'Energy = ', p % E + write(ou,*) ' x index = ' // int_to_str(p % index_x) + write(ou,*) ' y index = ' // int_to_str(p % index_y) + write(ou,*) ' IE = ' // int_to_str(p % IE) + write(ou,100) 'Interpolation factor = ', p % interp - if (p % surface > 0) then - s => surfaces(p % surface) - write(ou,*) ' Surface = ' // int_to_str(s % uid) - else - write(ou,*) ' Surface = None' - end if + if (p % cell > 0) then + c => cells(p % cell) + write(ou,*) ' Cell = ' // int_to_str(c % uid) + else + write(ou,*) ' Cell not determined' + end if - u => universes(p % universe) - write(ou,*) ' Universe = ' // int_to_str(u % uid) - write(ou,*) + if (p % surface > 0) then + s => surfaces(p % surface) + write(ou,*) ' Surface = ' // int_to_str(s % uid) + else + write(ou,*) ' Surface = None' + end if - ! Format for a single real - 100 format (5X,A,G10.3) + u => universes(p % universe) + write(ou,*) ' Universe = ' // int_to_str(u % uid) + write(ou,*) - nullify(c) - nullify(s) - nullify(u) + ! Format for a single real +100 format (5X,A,G10.3) - end subroutine print_particle + nullify(c) + nullify(s) + nullify(u) + + end subroutine print_particle !===================================================================== ! PRINT_REACTION displays the attributes of a reaction !===================================================================== - subroutine print_reaction(rxn) + subroutine print_reaction(rxn) - type(AceReaction), pointer :: rxn + type(AceReaction), pointer :: rxn - write(ou,*) 'Reaction ' // reaction_name(rxn % MT) - write(ou,*) ' MT = ' // int_to_str(rxn % MT) - write(ou,100) 'Q-value = ', rxn % Q_value - write(ou,*) ' TY = ' // int_to_str(rxn % TY) - write(ou,*) ' Starting index = ' // int_to_str(rxn % IE) - if (rxn % has_energy_dist) then - write(ou,*) ' Energy: Law ' // int_to_str(rxn % edist % law) - end if - write(ou,*) + write(ou,*) 'Reaction ' // reaction_name(rxn % MT) + write(ou,*) ' MT = ' // int_to_str(rxn % MT) + write(ou,100) 'Q-value = ', rxn % Q_value + write(ou,*) ' TY = ' // int_to_str(rxn % TY) + write(ou,*) ' Starting index = ' // int_to_str(rxn % IE) + if (rxn % has_energy_dist) then + write(ou,*) ' Energy: Law ' // int_to_str(rxn % edist % law) + end if + write(ou,*) - ! Format for a single real - 100 format (5X,A,G10.3) + ! Format for a single real +100 format (5X,A,G10.3) - end subroutine print_reaction + end subroutine print_reaction !===================================================================== ! PRINT_CELL displays the attributes of a cell !===================================================================== - subroutine print_cell(c) + subroutine print_cell(c) - type(Cell), pointer :: c + type(Cell), pointer :: c - integer :: temp - integer :: i - type(Universe), pointer :: u => null() - type(Material), pointer :: m => null() - character(250) :: string + integer :: temp + integer :: i + type(Universe), pointer :: u => null() + type(Lattice), pointer :: l => null() + type(Material), pointer :: m => null() + character(250) :: string - 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 - write(ou,*) ' Parent Cell = ' // int_to_str(c % parent) - string = "" - do i = 1, c % n_surfaces - select case (c % surfaces(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 % surfaces(i)) - end select - end do - write(ou,*) ' Surface Specification:' // trim(string) - write(ou,*) + 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) + select case (c % type) + case (CELL_NORMAL) + write(ou,*) ' Fill = NONE' + case (CELL_FILL) + u => universes(c % fill) + write(ou,*) ' Fill = Universe ' // int_to_str(u % uid) + case (CELL_LATTICE) + l => lattices(c % fill) + write(ou,*) ' Fill = Lattice ' // int_to_str(l % uid) + end select + if (c % material == 0) then + write(ou,*) ' Material = NONE' + else + m => materials(c % material) + write(ou,*) ' Material = ' // int_to_str(m % uid) + end if + write(ou,*) ' Parent Cell = ' // int_to_str(c % parent) + string = "" + do i = 1, c % n_surfaces + select case (c % surfaces(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 % surfaces(i)) + end select + end do + write(ou,*) ' Surface Specification:' // trim(string) + write(ou,*) - ! nullify associated pointers - nullify(u) - nullify(m) + ! nullify associated pointers + nullify(u) + nullify(m) - end subroutine print_cell + end subroutine print_cell !===================================================================== ! PRINT_UNIVERSE displays the attributes of a universe !===================================================================== - subroutine print_universe(univ) + subroutine print_universe(univ) - type(Universe), pointer :: univ + type(Universe), pointer :: univ - integer :: i - character(250) :: string - type(Cell), pointer :: c => null() + integer :: i + character(250) :: string + type(Cell), pointer :: c => null() - 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,*) + 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,*) - nullify(c) + nullify(c) - end subroutine print_universe + end subroutine print_universe + +!===================================================================== +! PRINT_LATTICE displays the attributes of a lattice +!===================================================================== + + subroutine print_lattice(lat) + + type(Lattice), pointer :: lat + + write(ou,*) 'Lattice ' // int_to_str(lat % uid) + write(ou,*) ' n_x = ' // int_to_str(lat % n_x) + write(ou,*) ' n_y = ' // int_to_str(lat % n_y) + write(ou,100) 'x0 = ', lat % x0 + write(ou,100) 'y0 = ', lat % y0 + write(ou,100) 'width_x = ', lat % width_x + write(ou,100) 'width_y = ', lat % width_y + write(ou,*) + + ! Format for a single real +100 format (5X,A,G10.3) + + end subroutine print_lattice !===================================================================== ! PRINT_SURFACE displays the attributes of a surface !===================================================================== - subroutine print_surface(surf) + subroutine print_surface(surf) - type(Surface), pointer :: surf + type(Surface), pointer :: surf - integer :: i - character(80) :: string + integer :: i + character(80) :: string - write(ou,*) 'Surface ' // int_to_str(surf % uid) - select case (surf % type) - case (SURF_PX) - string = "X Plane" - case (SURF_PY) - string = "Y Plane" - case (SURF_PZ) - string = "Z Plane" - case (SURF_PLANE) - string = "Plane" - case (SURF_CYL_X) - string = "X Cylinder" - case (SURF_CYL_Y) - string = "Y Cylinder" - case (SURF_CYL_Z) - string = "Z Cylinder" - case (SURF_SPHERE) - string = "Sphere" - case (SURF_BOX_X) - case (SURF_BOX_Y) - case (SURF_BOX_Z) - case (SURF_BOX) - case (SURF_GQ) - string = "General Quadratic" - end select - write(ou,*) ' Type = ' // trim(string) - write(ou,*) ' Coefficients = ', surf % coeffs + write(ou,*) 'Surface ' // int_to_str(surf % uid) + select case (surf % type) + case (SURF_PX) + string = "X Plane" + case (SURF_PY) + string = "Y Plane" + case (SURF_PZ) + string = "Z Plane" + case (SURF_PLANE) + string = "Plane" + case (SURF_CYL_X) + string = "X Cylinder" + case (SURF_CYL_Y) + string = "Y Cylinder" + case (SURF_CYL_Z) + string = "Z Cylinder" + case (SURF_SPHERE) + string = "Sphere" + case (SURF_BOX_X) + case (SURF_BOX_Y) + case (SURF_BOX_Z) + case (SURF_BOX) + case (SURF_GQ) + string = "General Quadratic" + end select + write(ou,*) ' Type = ' // trim(string) + write(ou,*) ' Coefficients = ', surf % coeffs - string = "" - if (allocated(surf % neighbor_pos)) then - do i = 1, size(surf % neighbor_pos) - string = trim(string) // ' ' // int_to_str(surf % neighbor_pos(i)) - end do - end if - write(ou,*) ' Positive Neighbors = ', trim(string) + string = "" + if (allocated(surf % neighbor_pos)) then + do i = 1, size(surf % neighbor_pos) + string = trim(string) // ' ' // int_to_str(surf % neighbor_pos(i)) + end do + end if + write(ou,*) ' Positive Neighbors = ', trim(string) - string = "" - if (allocated(surf % neighbor_neg)) then - do i = 1, size(surf % neighbor_neg) - string = trim(string) // ' ' // int_to_str(surf % neighbor_neg(i)) - end do - end if - write(ou,*) ' Negative Neighbors =', trim(string) + string = "" + if (allocated(surf % neighbor_neg)) then + do i = 1, size(surf % neighbor_neg) + string = trim(string) // ' ' // int_to_str(surf % neighbor_neg(i)) + end do + end if + write(ou,*) ' Negative Neighbors =', trim(string) + write(ou,*) - end subroutine print_surface + end subroutine print_surface !===================================================================== ! PRINT_MATERIAL displays the attributes of a material !===================================================================== - subroutine print_material(mat) + subroutine print_material(mat) - type(Material), pointer :: mat + type(Material), pointer :: mat - integer :: i - integer :: n_lines - type(AceContinuous), pointer :: table - character(250) :: string + integer :: i + integer :: n_lines + type(AceContinuous), pointer :: table + character(250) :: string - write(ou,*) 'Material ' // int_to_str(mat % uid) - ! Make string of all isotopes - string = "" - do i = 1, mat % n_isotopes - table => xs_continuous(mat % table(i)) - string = trim(string) // ' ' // table % name - end do - ! Print isotopes with word wrap - ! TODO: Change this to generic word wrap subroutine? - n_lines = (len_trim(string)-1)/75 + 1 - do i = 1, n_lines - if ( i == 1 ) then - write(ou, fmt='(5X,A75)') 'Isotopes =' // string(75*(i-1)+1:75*i) - else - write(ou, fmt='(5X,A75)') string(75*(i-1)+1:75*i) - end if - end do - write(ou,'(5X,A,G12.4,A)') 'Atom Density = ', mat % atom_density, & - & ' atom/b-cm' - write(ou,*) + write(ou,*) 'Material ' // int_to_str(mat % uid) + ! Make string of all isotopes + string = "" + do i = 1, mat % n_isotopes + table => xs_continuous(mat % table(i)) + string = trim(string) // ' ' // table % name + end do + ! Print isotopes with word wrap + ! TODO: Change this to generic word wrap subroutine? + n_lines = (len_trim(string)-1)/75 + 1 + do i = 1, n_lines + if ( i == 1 ) then + write(ou, fmt='(5X,A75)') 'Isotopes =' // string(75*(i-1)+1:75*i) + else + write(ou, fmt='(5X,A75)') string(75*(i-1)+1:75*i) + end if + end do + write(ou,'(5X,A,G12.4,A)') 'Atom Density = ', mat % atom_density, & + & ' atom/b-cm' + write(ou,*) - nullify(table) + nullify(table) - end subroutine print_material + end subroutine print_material !===================================================================== ! PRINT_SUMMARY displays the attributes of all cells, universes, @@ -488,59 +519,73 @@ module output ! debugging! !===================================================================== - subroutine print_summary() + subroutine print_summary() - type(Surface), pointer :: s - type(Cell), pointer :: c - type(Universe), pointer :: u - type(Material), pointer :: m - integer :: i + type(Surface), pointer :: s => null() + type(Cell), pointer :: c => null() + type(Universe), pointer :: u => null() + type(Lattice), pointer :: l => null() + type(Material), pointer :: m => null() + integer :: i - ! 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 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 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 - write(ou,*) '=============================================' - write(ou,*) '=> SURFACE SUMMARY <=' - write(ou,*) '=============================================' - write(ou,*) - do i = 1, n_surfaces - s => surfaces(i) - call print_surface(s) - end do + ! print summary of lattices + if (n_lattices > 0) then + write(ou,*) '=============================================' + write(ou,*) '=> LATTICE SUMMARY <=' + write(ou,*) '=============================================' + write(ou,*) + do i = 1, n_lattices + l => lattices(i) + call print_lattice(l) + end do + end if - ! print summary of materials - write(ou,*) '=============================================' - write(ou,*) '=> MATERIAL SUMMARY <=' - write(ou,*) '=============================================' - write(ou,*) - do i = 1, n_materials - m => materials(i) - call print_material(m) - end do + ! print summary of surfaces + write(ou,*) '=============================================' + write(ou,*) '=> SURFACE SUMMARY <=' + write(ou,*) '=============================================' + write(ou,*) + do i = 1, n_surfaces + s => surfaces(i) + call print_surface(s) + end do - nullify(s) - nullify(c) - nullify(u) - nullify(m) + ! print summary of materials + write(ou,*) '=============================================' + write(ou,*) '=> MATERIAL SUMMARY <=' + write(ou,*) '=============================================' + write(ou,*) + do i = 1, n_materials + m => materials(i) + call print_material(m) + end do - end subroutine print_summary + nullify(s) + nullify(c) + nullify(u) + nullify(l) + nullify(m) + + end subroutine print_summary end module output diff --git a/src/physics.f90 b/src/physics.f90 index 385915a6e..dd63b0694 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -1,7 +1,8 @@ module physics use global - use geometry, only: find_cell, dist_to_boundary, cross_boundary + use geometry, only: find_cell, dist_to_boundary, cross_surface, & + & cross_lattice use types, only: Particle use mcnp_random, only: rang use output, only: error, warning, message, print_particle @@ -29,6 +30,7 @@ contains real(8) :: Sigma ! total cross-section real(8) :: f ! interpolation factor logical :: found_cell ! found cell which particle is in? + logical :: in_lattice ! is surface crossing in lattice? character(250) :: msg ! output/error message type(Universe), pointer :: univ @@ -67,7 +69,7 @@ contains Sigma = (1-f)*cMaterial%total_xs(IE) + f*cMaterial%total_xs(IE+1) ! Find the distance to the nearest boundary - call dist_to_boundary(p, d_to_boundary, surf) + call dist_to_boundary(p, d_to_boundary, surf, in_lattice) ! Sample a distance to collision d_to_collision = -log(rang()) / Sigma @@ -77,13 +79,19 @@ contains ! Advance particle 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 % surface = surf p % cell = 0 - call cross_boundary(p) + if (in_lattice) then + p % surface = 0 + call cross_lattice(p) + else + p % surface = surf + call cross_surface(p) + end if else ! collision p % surface = 0 diff --git a/src/source.f90 b/src/source.f90 index 80ddbcdd4..6c64eaf2a 100644 --- a/src/source.f90 +++ b/src/source.f90 @@ -65,6 +65,7 @@ contains r = (/ (rang(), k = 1,3) /) p % uid = j p % xyz = p_min + r*(p_max - p_min) + p % xyz_local = p % xyz ! angle phi = TWO*PI*rang() @@ -78,8 +79,11 @@ contains p % cell = 0 p % surface = 0 p % universe = 0 + p % lattice = 0 p % wgt = ONE p % alive = .true. + p % index_x = 0 + p % index_y = 0 ! sample energy from Watt fission energy spectrum for U-235 do @@ -160,10 +164,11 @@ contains do i = 1, n_sites j = index + i - 1 - source_bank(j) % uid = temp_bank(i) % uid - source_bank(j) % xyz = temp_bank(i) % xyz - source_bank(j) % uvw = temp_bank(i) % uvw - source_bank(j) % E = temp_bank(i) % E + source_bank(j) % uid = temp_bank(i) % uid + source_bank(j) % xyz = temp_bank(i) % xyz + source_bank(j) % xyz_local = temp_bank(i) % xyz + source_bank(j) % uvw = temp_bank(i) % uvw + source_bank(j) % E = temp_bank(i) % E ! set defaults source_bank(j) % type = NEUTRON diff --git a/src/types.f90 b/src/types.f90 index decc03908..edf5a307c 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -3,7 +3,7 @@ module types implicit none !===================================================================== -! UNIVERSE +! UNIVERSE defines a geometry that fills all phase space !===================================================================== type Universe @@ -18,14 +18,21 @@ module types end type Universe !===================================================================== -! LATTICE +! LATTICE is an ordered array of elements (either rectangular, +! hexagonal, or triangular) !===================================================================== type Lattice - integer :: uid ! Unique ID - integer :: type ! Type of lattice (square, hex, etc) - integer :: level ! Level of lattice - real(8) :: pitch ! Lattice pitch in cm + integer :: uid ! Universe number for lattice + integer :: type ! Type of lattice (rectangular, hex, etc) + integer :: level ! Level of lattice + integer :: n_x ! number of lattice cells in x-direction + integer :: n_y ! number of lattice cells in y-direction + real(8) :: x0 ! x-coordinate of lattice origin + real(8) :: y0 ! y-coordinate of lattice origin + real(8) :: width_x ! width of lattice cell + real(8) :: width_y ! width of lattice cell + integer, allocatable :: element(:,:) ! specified universes end type Lattice !===================================================================== @@ -67,19 +74,23 @@ module types !===================================================================== 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) :: wgt ! particle weight - 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 - logical :: alive ! is particle alive? - integer :: n_coll ! # of collisions + integer :: uid ! Unique ID + integer :: type ! Particle type (n, p, e, etc) + real(8) :: xyz(3) ! location + real(8) :: xyz_local(3) ! local location (after transformations) + real(8) :: uvw(3) ! directional cosines + real(8) :: wgt ! particle weight + 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 :: lattice ! current lattice + integer :: surface ! current surface + integer :: index_x ! lattice index for x direction + integer :: index_y ! lattice index for y direction + logical :: alive ! is particle alive? + integer :: n_coll ! # of collisions end type Particle !=====================================================================