From 4d8cca60ece70da34e23e3b20ca6b79bb8fac8d8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 4 Mar 2011 17:26:26 +0000 Subject: [PATCH] Universe capability mostly complete now. Initial lattice routines and structures added. --- ChangeLog | 12 +++++ src/fileio.f90 | 64 ++++++++++++++++++++++++-- src/geometry.f90 | 70 +++++++++++++++++++++------- src/global.f90 | 4 ++ src/input_sample | 6 +-- src/main.f90 | 10 +++- src/output.f90 | 116 +++++++++++++++++++++++++++++++++++++++++++++++ src/physics.f90 | 56 ++++++++++++++++++----- src/types.f90 | 16 ++++++- 9 files changed, 315 insertions(+), 39 deletions(-) diff --git a/ChangeLog b/ChangeLog index a9fd23cf0..41f0a22ec 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,15 @@ +2011-03-04 Paul Romano + + * fileio.f90: Added build_universe routine which sets pointers + betweens cells and parent cells and determines level of each + universe. Added stub for read_lattice. + * geometry.f90: Changed cross_boundary and dist_to_boundary to + take into account universe boundaries as well as cell boundaries. + * global.f90: Added initial variables for lattices. + * output.f90: Added print_surface and print_material. + * physics.f90: Added find_energy_index routine. + * types.f90: Initial Lattice type. + 2011-02-27 Paul Romano * ace.f90: Fixed minor bugs associated with reading angular diff --git a/src/fileio.f90 b/src/fileio.f90 index 4c1cf8148..7430af6c8 100644 --- a/src/fileio.f90 +++ b/src/fileio.f90 @@ -214,6 +214,7 @@ contains integer :: index_cell integer :: index_universe integer :: index_surface + integer :: index_lattice integer :: index_material integer :: index_source @@ -223,6 +224,7 @@ contains index_cell = 0 index_universe = 0 index_surface = 0 + index_lattice = 0 index_material = 0 index_source = 0 @@ -251,6 +253,11 @@ contains index_surface = index_surface + 1 call read_surface(index_surface, words, n) + case ('lattice') + ! Read data from lattice entry + index_lattice = index_lattice + 1 + call read_lattice(index_lattice, words, n) + case ('material') index_material = index_material + 1 call read_material(index_material, words, n) @@ -288,7 +295,7 @@ contains subroutine adjust_indices() - type(Cell), pointer :: c + type(Cell), pointer :: c => null() integer :: i, j integer :: index @@ -328,6 +335,41 @@ contains end subroutine adjust_indices +!===================================================================== +! BUILD_UNIVERSE determines what level each universe is at and +! determines what the parent cell of each cell in a subuniverse is. +!===================================================================== + + recursive subroutine build_universe(univ, parent, level) + + type(Universe), pointer :: univ + integer, intent(in) :: parent + integer, intent(in) :: level + + integer :: i + integer :: index + type(Cell), pointer :: c => null() + type(Universe), pointer :: subuniverse => null() + + ! set level of the universe + univ % level = level + + ! loop over all cells in the universe + do i = 1, univ % n_cells + index = univ % cells(i) + c => cells(index) + c%parent = parent + + ! if this cell is filled with another universe, recursive call + ! this subroutine + if (c % fill > 0) then + subuniverse => universes(c % fill) + call build_universe(subuniverse, index, level + 1) + end if + end do + + end subroutine build_universe + !===================================================================== ! READ_CELL parses the data on a cell card. !===================================================================== @@ -508,7 +550,19 @@ contains end subroutine read_surface !===================================================================== -! READ_SOURCE parses the data on a source card. +! READ_LATTICE parses the data on a lattice entry. +!===================================================================== + + subroutine read_lattice(index, words, n_words) + + integer, intent(in) :: index + character(*), intent(in) :: words(max_words) + integer, intent(in) :: n_words + + end subroutine read_lattice + +!===================================================================== +! READ_SOURCE parses the data on a source entry. !===================================================================== subroutine read_source(words, n_words) @@ -563,7 +617,7 @@ contains integer :: i integer :: ioError integer :: n_isotopes - type(Material), pointer :: mat + type(Material), pointer :: mat => null() if (mod(n_words,2) == 0 .or. n_words < 5) then msg = "Invalid number of arguments for material: " // words(2) @@ -606,8 +660,8 @@ contains subroutine normalize_ao() - type(xsData), pointer :: iso - type(Material), pointer :: mat + type(xsData), pointer :: iso => null() + type(Material), pointer :: mat => null() integer :: index integer :: i, j real(8) :: sum_percent diff --git a/src/geometry.f90 b/src/geometry.f90 index 1be471cd1..0bb71b092 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -142,10 +142,12 @@ contains integer :: i integer :: index_cell character(250) :: msg + logical :: found + type(Universe), pointer :: lower_univ => null() surf => surfaces(abs(neut%surface)) if (verbosity >= 10) then - msg = " Crossing surface " // trim(int_to_str(neut%surface)) + msg = " Crossing surface " // trim(int_to_str(surf%uid)) call message(msg, 10) end if @@ -159,29 +161,49 @@ contains return end if - if (neut%surface < 0 .and. allocated(surf%neighbor_pos)) then + if (neut%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 - neut%cell = index_cell - cCell => cells(index_cell) - cMaterial => materials(cCell%material) + if (c % fill > 0) then + lower_univ => universes(c % fill) + call find_cell(lower_univ, neut, found) + if (.not. found) then + msg = "Could not locate neutron in universe: " + call error(msg) + end if + else + ! set current pointers + neut%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 (neut%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 - neut%cell = index_cell - cCell => cells(index_cell) - cMaterial => materials(cCell%material) + if (c % fill > 0) then + lower_univ => universes(c % fill) + call find_cell(lower_univ, neut, found) + if (.not. found) then + msg = "Could not locate neutron in universe: " + call error(msg) + end if + else + ! set current pointers + neut%cell = index_cell + cCell => c + cMaterial => materials(cCell%material) + end if return end if end do @@ -193,7 +215,7 @@ contains c => cells(i) if (cell_contains(c, neut)) then neut%cell = i - cCell => cells(i) + cCell => c cMaterial => materials(cCell%material) return end if @@ -209,7 +231,8 @@ contains !===================================================================== ! DIST_TO_BOUNDARY calculates the distance to the nearest boundary of ! the cell 'cl' for a particle 'neut' traveling in a certain -! direction. +! 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) @@ -220,9 +243,10 @@ contains integer, optional, intent(in) :: other_cell type(Cell), pointer :: cell_p => null() + type(Cell), pointer :: parent_p => null() type(Surface), pointer :: surf_p => null() integer :: i - integer :: n_boundaries + integer :: n_boundaries, n1, n2 integer, allocatable :: expression(:) integer :: index_surf integer :: current_surf @@ -242,10 +266,24 @@ contains current_surf = neut%surface - dist = INFINITY - n_boundaries = size(cell_p%boundary_list) + ! determine number of surfaces to check + n1 = cell_p % n_items + n2 = 0 + if (cell_p % parent > 0) then + parent_p => cells(cell_p % parent) + n2 = parent_p % n_items + end if + n_boundaries = n1 + n2 + + ! allocate space and assign expression allocate(expression(n_boundaries)) - expression = cell_p%boundary_list + expression(1:n1) = cell_p%boundary_list + if (cell_p % parent > 0) then + expression(n1+1:n1+n2) = parent_p % boundary_list + end if + + ! loop over all surfaces + dist = INFINITY do i = 1, n_boundaries x = neut%xyz(1) y = neut%xyz(2) @@ -254,7 +292,7 @@ contains v = neut%uvw(2) w = neut%uvw(3) - ! check for coincident surfaec + ! check for coincident surface index_surf = expression(i) if (index_surf == current_surf) cycle diff --git a/src/global.f90 b/src/global.f90 index 0a9c15c93..d9231ffe0 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -7,18 +7,21 @@ module global ! Main arrays for cells, surfaces, materials type(Cell), allocatable, target :: cells(:) type(Universe), allocatable, target :: universes(:) + type(Lattice), allocatable, target :: lattices(:) 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_lattices ! # of lattices 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 :: lattice_dict type(DictionaryII), pointer :: surface_dict type(DictionaryII), pointer :: material_dict type(DictionaryII), pointer :: isotope_dict @@ -34,6 +37,7 @@ module global ! Current cell, surface, material type(Cell), pointer :: cCell type(Universe), pointer :: cUniverse + type(Lattice), pointer :: cLattice type(Surface), pointer :: cSurface type(Material), pointer :: cMaterial diff --git a/src/input_sample b/src/input_sample index e067a1dd8..b4319adf7 100644 --- a/src/input_sample +++ b/src/input_sample @@ -5,9 +5,9 @@ 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 +surface 1 cylz 0 0 3 +surface 2 cylz 0 0 5 +surface 3 cylz 0 0 8 material 40 -4.5 & 92238.03c 1.0 diff --git a/src/main.f90 b/src/main.f90 index 771828ae5..19e966ff2 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -2,7 +2,7 @@ program main use global use fileio, only: read_input, read_command_line, read_count, & - & normalize_ao + & normalize_ao, build_universe use output, only: title, echo_input, message, warning, error, & & print_summary use geometry, only: sense, cell_contains, neighbor_lists @@ -17,6 +17,7 @@ program main character(16) :: filename character(250) :: msg + type(Universe), pointer :: univ ! Print the OpenMC title and version/date/time information call title() @@ -37,6 +38,11 @@ program main call read_count(path_input) call read_input(path_input) + ! determine at which level universes are and link cells to parenting + ! cells + univ => universes(BASE_UNIVERSE) + call build_universe(univ, 0, 0) + ! After reading input and basic geometry setup is complete, build ! lists of neighboring cells for efficient tracking call neighbor_lists() @@ -66,7 +72,7 @@ program main call init_source() ! start problem - surfaces(2)%bc = BC_VACUUM + surfaces(3)%bc = BC_VACUUM call run_problem() ! deallocate arrays diff --git a/src/output.f90 b/src/output.f90 index ad74f5480..dfc2f5b28 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -254,6 +254,7 @@ module output 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_items select case (c % boundary_List(i)) @@ -305,6 +306,105 @@ module output end subroutine print_universe +!===================================================================== +! PRINT_SURFACE displays the attributes of a surface +!===================================================================== + + subroutine print_surface(surf) + + type(Surface), pointer :: surf + + integer :: ou + integer :: i + character(80) :: string + + ou = OUTPUT_UNIT + + 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 = ', 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 =', string + + end subroutine print_surface + +!===================================================================== +! PRINT_MATERIAL displays the attributes of a material +!===================================================================== + + subroutine print_material(mat) + + type(Material), pointer :: mat + + integer :: ou + integer :: i + integer :: n_lines + type(AceContinuous), pointer :: table + character(250) :: string + + ou = OUTPUT_UNIT + + 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,*) + + end subroutine print_material + !===================================================================== ! PRINT_SUMMARY displays the attributes of all cells, universes, ! surfaces and materials read in the input file. Very useful for @@ -343,8 +443,24 @@ module output 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 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 diff --git a/src/physics.f90 b/src/physics.f90 index 6786dd982..d9035de62 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, print_universe + use output, only: error, message use search, only: binary_search use endf, only: reaction_name @@ -38,7 +38,6 @@ contains if (neut%cell == 0) then 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 @@ -48,8 +47,6 @@ contains end if end if - return - if (verbosity >= 10) then msg = "=== Particle " // trim(int_to_str(neut%uid)) // " ===" call message(msg, 10) @@ -63,13 +60,12 @@ contains neut%E = watt_spectrum(0.988_8, 2.249_8) ! find energy index, interpolation factor - IE = binary_search(e_grid, n_grid, neut%E) - f = (neut%E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE)) + call find_energy_index(neut) + IE = neut % IE + f = neut % interp ! Determine material total cross-section Sigma = f*cMaterial%total_xs(IE) + (1-f)*cMaterial%total_xs(IE+1) - neut%IE = IE - neut%interp = f do while (neut%alive) @@ -78,6 +74,8 @@ contains ! Sample a distance to collision d_to_collision = -log(rang()) / 1.0 ! Sigma + + ! Select smaller of the two distances distance = min(d_to_boundary, d_to_collision) ! Advance neutron @@ -98,6 +96,44 @@ contains end subroutine transport +!===================================================================== +! FIND_ENERGY_INDEX determines the index on the union energy grid and +! the interpolation factor for a particle at a certain energy +!===================================================================== + + subroutine find_energy_index(neut) + + type(Neutron), pointer :: neut + + integer :: IE + real(8) :: E + real(8) :: interp + + ! copy particle's energy + E = neut%E + + ! if particle's energy is outside of energy grid range, set to + ! first or last index. Otherwise, do a binary search through the + ! union energy grid. + if (E < e_grid(1)) then + IE = 1 + elseif (E > e_grid(n_grid)) then + IE = n_grid + else + IE = binary_search(e_grid, n_grid, E) + end if + + ! calculate the interpolation factor -- note this will be outside + ! of [0,1) for a particle outside the energy range of the union + ! grid + interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE)) + + ! set particle attributes + neut % IE = IE + neut % interp = interp + + end subroutine find_energy_index + !===================================================================== ! COLLISION !===================================================================== @@ -247,9 +283,7 @@ contains neut%uvw(3) = vz/vel ! find energy index, interpolation factor - IE = binary_search(e_grid, n_grid, E) - neut%IE = IE - neut%interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE)) + call find_energy_index(neut) end subroutine elastic_scatter diff --git a/src/types.f90 b/src/types.f90 index ceaa2ad6a..6f83b5add 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -17,6 +17,17 @@ module types real(8) :: z0 end type Universe +!===================================================================== +! LATTICE +!===================================================================== + + type Lattice + integer :: uid + integer :: type + integer :: level + real(8) :: pitch + end type Lattice + !===================================================================== ! SURFACE type defines a first- or second-order surface that can be ! used to construct closed volumes (cells) @@ -38,8 +49,9 @@ module types type Cell integer :: uid integer :: type - integer :: universe - integer :: fill + 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(:)