crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Universe capability mostly complete now. Initial lattice routines and structures added.

Browse source
This commit is contained in:
Paul Romano 2011-03-04 17:26:26 +00:00
parent 03a2634fb4
commit 4d8cca60ec
9 changed files with 315 additions and 39 deletions
Download patch file Download diff file Expand all files Collapse all files

12
ChangeLog
View file

@ -1,3 +1,15 @@
2011-03-04 Paul Romano <romano7@mit.edu>
* 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 <romano7@mit.edu>
* ace.f90: Fixed minor bugs associated with reading angular

64
src/fileio.f90
View file

@ -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

70
src/geometry.f90
View file

@ -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

4
src/global.f90
View file

@ -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

6
src/input_sample
View file

@ -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

10
src/main.f90
View file

@ -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

116
src/output.f90
View file

@ -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

56
src/physics.f90
View file

@ -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

16
src/types.f90
View file

@ -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(:)