diff --git a/ChangeLog b/ChangeLog index b4656841b..a58ede0fd 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,16 @@ +2011-01-23 Paul Romano + + * types.f90: Added cell, surface, and alive attributes to Neutron + * fileio.f90: External source no longer array + * geometry.f90: Changed cell_contains slightly so that it can + evaluate simple cell logical expressions. Added find_cell and + cross_boundary subroutines. Changed dist_to_boundary to also give + back closest surface + * global.f90: Added set_defaults subroutine + * source.f90: Isotropic angle is now sampled in init_source, fixed + sampling in box, checks for external source type. Added + get_source_particle function. + 2011-01-22 Paul Romano * input_sample: Added a sample input file diff --git a/src/Makefile b/src/Makefile index 2e11fd2a0..1649c82a7 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,7 +1,8 @@ program = openmc src = main.f90 types.f90 global.f90 fileio.f90 output.f90 \ - string.f90 geometry.f90 mcnp_random.f90 source.f90 + string.f90 geometry.f90 mcnp_random.f90 source.f90 \ + physics.f90 objects = $(src:.f90=.o) #-------------------------------------------------------------------- @@ -34,7 +35,8 @@ global.o: types.o string.o: global.o output.o output.o: global.o source.o: global.o mcnp_random.o +physics.o: types.o global.o geometry.o: types.o global.o output.o string.o fileio.o: types.o global.o string.o output.o main.o: global.o fileio.o output.o geometry.o mcnp_random.o \ - source.o + source.o physics.o diff --git a/src/fileio.f90 b/src/fileio.f90 index 0bb37f984..d343c8989 100644 --- a/src/fileio.f90 +++ b/src/fileio.f90 @@ -57,10 +57,9 @@ contains integer :: readError integer :: n - n_cells = 0 - n_surfaces = 0 - n_materials = 0 - n_sources = 0 + n_cells = 0 + n_surfaces = 0 + n_materials = 0 open(file=filename, unit=in, status='old', action='read') @@ -77,13 +76,11 @@ contains n_surfaces = n_surfaces + 1 case ( 'material' ) n_materials = n_materials + 1 - case ( 'source' ) - n_sources = n_sources + 1 end select end do - ! Check to make sure there are cells, surface, materials, and - ! sources defined for the problem + ! Check to make sure there are cells, surface, and materials + ! defined for the problem if ( n_cells == 0 ) then msg = "No cells specified!" call error( msg ) @@ -93,9 +90,6 @@ contains elseif ( n_materials == 0 ) then msg = "No materials specified!" call error( msg ) - elseif ( n_sources == 0 ) then - msg = "No source specified!" - call error( msg ) end if close(unit=in) @@ -104,7 +98,6 @@ contains allocate( cells(n_cells) ) allocate( surfaces(n_surfaces) ) allocate( materials(n_materials) ) - allocate( sources(n_sources) ) end subroutine read_count @@ -157,8 +150,7 @@ contains case ( 'material' ) case ( 'source' ) - index_source = index_source + 1 - call read_source( index_source, words, n ) + call read_source( words, n ) case ( 'criticality' ) call read_criticality( words, n ) @@ -312,50 +304,39 @@ contains ! READ_SOURCE parses the data on a source card. !===================================================================== - subroutine read_source( index, words, n_words ) + subroutine read_source( words, n_words ) - integer, intent(in) :: index character(*), intent(in) :: words(max_words) integer, intent(in) :: n_words character(250) :: msg character(32) :: word - type(ExtSource), pointer :: this_source => null() integer :: readError integer :: values_reqd integer :: i - this_source => sources(index) - - ! Read source identifier - read(words(2), fmt='(I8)', iostat=readError) this_source%uid - if ( readError > 0 ) then - msg = "Invalid surface name: " // words(2) - call error( msg ) - end if - ! Read source type - word = words(3) + word = words(2) call lower_case(word) select case ( trim(word) ) case ( 'box' ) - this_source%type = SRC_BOX + external_source%type = SRC_BOX values_reqd = 6 case default - msg = "Invalid source type: " // words(3) + msg = "Invalid source type: " // words(2) call error( msg ) end select ! Make sure there are enough values for this source type - if ( n_words-3 < values_reqd ) then - msg = "Not enough values for source: " // words(2) + if ( n_words-2 < values_reqd ) then + msg = "Not enough values for source of type: " // words(2) call error( msg ) end if ! Read list of surfaces - allocate( this_source%values(n_words-3) ) - do i = 1, n_words-3 - this_source%values(i) = str_to_real(words(i+3)) + allocate( external_source%values(n_words-2) ) + do i = 1, n_words-2 + external_source%values(i) = str_to_real(words(i+2)) end do end subroutine read_source diff --git a/src/geometry.f90 b/src/geometry.f90 index 465a8982e..5f41246e2 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -12,67 +12,168 @@ contains ! CELL_CONTAINS determines whether a given point is inside a cell !===================================================================== - subroutine cell_contains( c, xyz, on_surface ) + function cell_contains( c, xyz ) type(Cell), intent(in) :: c real(8), intent(in) :: xyz(3) - logical, intent(in) :: on_surface + logical :: cell_contains integer, allocatable :: expression(:) integer :: specified_sense integer :: actual_sense - integer :: n + integer :: n_boundaries integer :: i, j integer :: surf_num type(Surface), pointer :: surf => null() logical :: surface_found character(250) :: msg - if ( on_surface ) then - n = size(c%boundary_list) - allocate( expression(n) ) - expression = c%boundary_list - do i = 1,n - ! Don't change logical operator - if ( expression(i) >= OP_DIFFERENCE ) then - cycle - end if + n_boundaries = size(c%boundary_list) + allocate( expression(n_boundaries) ) + expression = c%boundary_list + do i = 1, n_boundaries + surface_found = .false. - ! Lookup surface - surf_num = abs(expression(i)) - ! TODO: replace this loop with a hash since this lookup is O(N) - do j = 1, n_surfaces - surf => surfaces(j) - if ( surf%uid == surf_num ) then - surface_found = .true. - exit - end if - end do + ! Don't change logical operator + if ( expression(i) >= OP_DIFFERENCE ) then + cycle + end if - ! Report error if can't find specified surface - if ( .not. surface_found ) then - deallocate( expression ) - msg = "Count not find surface " // trim(int_to_str(surf_num)) // & - & " specified on cell " // int_to_str(c%uid) - call error( msg ) + ! Lookup surface + surf_num = abs(expression(i)) + ! TODO: replace this loop with a hash since this lookup is O(N) + do j = 1, n_surfaces + surf => surfaces(j) + if ( surf%uid == surf_num ) then + surface_found = .true. + exit end if - - ! Compare sense of point to specified sense - specified_sense = sign(1,expression(i)) - call sense( surf, xyz, actual_sense ) - if ( actual_sense == specified_sense ) then - expression(i) = 1 - else - expression(i) = 0 + end do + + ! Report error if can't find specified surface + if ( .not. surface_found ) then + deallocate( expression ) + msg = "Count not find surface " // trim(int_to_str(surf_num)) // & + & " specified on cell " // int_to_str(c%uid) + call error( msg ) + end if + + ! Compare sense of point to specified sense + specified_sense = sign(1,expression(i)) + call sense( surf, xyz, actual_sense ) + if ( actual_sense == specified_sense ) then + expression(i) = 1 + else + expression(i) = 0 + end if + + end do + + ! TODO: Need to replace this with a 'lgeval' like subroutine that + ! can actually test expressions with unions and parentheses + if ( all( expression == 1 ) ) then + cell_contains = .true. + else + cell_contains = .false. + end if + + ! Free up memory from expression + deallocate( expression ) + + end function cell_contains + +!===================================================================== +! FIND_CELL determines what cell a source neutron is in +!===================================================================== + + subroutine find_cell( neut ) + + type(Neutron), pointer, intent(inout) :: neut + + logical :: found_cell + character(250) :: msg + integer :: i + + ! determine what region in + do i = 1, n_cells + if ( cell_contains(cells(i), neut%xyz) ) then + neut%cell = i + found_cell = .true. + exit + end if + end do + + ! if neutron couldn't be located, print error + if ( .not. found_cell ) then + write(msg, 100) "Could not locate cell for neutron at: ", neut%xyz +100 format (A,3ES11.3) + call error( msg ) + end if + + end subroutine find_cell + +!===================================================================== +! CROSS_BOUNDARY moves a neutron into a new cell +!===================================================================== + + subroutine cross_boundary( neut ) + + type(Neutron), pointer, intent(in) :: neut + + type(Surface), pointer :: surf + type(Cell), pointer :: c + integer :: i + integer :: index_cell + character(100) :: msg + + surf => surfaces(abs(neut%surface)) + + ! check for leakage + if ( surf%bc == BC_VACUUM ) then + neut%alive = .false. + return + end if + + if ( neut%surface < 0 ) 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%xyz) ) then + neut%cell = index_cell + return + end if + end do + elseif ( neut%surface > 0 ) 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%xyz) ) then + neut%cell = index_cell + return end if - end do end if - ! Need to deallocate expression - print *, expression + ! Couldn't find particle in neighboring cells, search through all + ! cells + do i = 1, size(cells) + c = cells(i) + if ( cell_contains(c, neut%xyz) ) then + neut%cell = i + return + end if + end do - end subroutine cell_contains + ! Couldn't find next cell anywhere! + msg = "After neutron crossed surface " // int_to_str(neut%surface) // & + & ", it could not be located in any cell and it did not leak." + call error( msg ) + + end subroutine cross_boundary !===================================================================== ! DIST_TO_BOUNDARY calculates the distance to the nearest boundary of @@ -80,17 +181,19 @@ contains ! direction. !===================================================================== - subroutine dist_to_boundary( cl, neut, dist ) + subroutine dist_to_boundary( neut, dist, surf, other_cell ) - type(Cell), intent(in) :: cl - type(Neutron), intent(in) :: neut - real(8), intent(out) :: dist + type(Neutron), intent(in) :: neut + real(8), intent(out) :: dist + integer, intent(out) :: surf + integer, optional, intent(in) :: other_cell + type(Cell), pointer :: cl => null() + type(Surface), pointer :: surf => null() integer :: i - integer :: n_cell_surfaces + integer :: n_boundaries integer, allocatable :: expression(:) integer :: surf_num - type(Surface), pointer :: surf => null() real(8) :: x,y,z,u,v,w real(8) :: d real(8) :: x0,y0,z0,r @@ -99,6 +202,12 @@ contains real(8) :: quad character(250) :: msg + if ( present(other_cell) ) then + cl => cells(other_cell) + else + cl => cells(neut%cell) + end if + x = neut%xyz(1) y = neut%xyz(2) z = neut%xyz(3) @@ -107,8 +216,8 @@ contains z = neut%uvw(3) dist = INFINITY - n_cell_surfaces = size(cl%boundary_list) - do i = 1, n_cell_surfaces + n_boundaries = size(cl%boundary_list) + do i = 1, n_boundaries surf_num = abs(expression(i)) if ( surf_num >= OP_DIFFERENCE ) cycle @@ -309,7 +418,10 @@ contains end select ! Check is calculated distance is new minimum - dist = min(d,dist) + if ( d < dist ) then + dist = d + surf = expression(i) + end if end do @@ -449,5 +561,3 @@ contains end subroutine sense end module geometry - - diff --git a/src/global.f90 b/src/global.f90 index 9e59178c4..b17ebfa8d 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -8,17 +8,18 @@ module global type(Cell), allocatable, target :: cells(:) type(Surface), allocatable, target :: surfaces(:) type(Material), allocatable, target :: materials(:) - type(ExtSource), allocatable, target :: sources(:) integer :: n_cells ! # of cells integer :: n_surfaces ! # of surfaces integer :: n_materials ! # of materials - integer :: n_sources ! # of sources ! Histories/cycles/etc for both external source and criticality integer :: n_particles ! # of particles (per cycle for criticality) integer :: n_cycles ! # of cycles integer :: n_inactive ! # of inactive cycles + ! External source + type(ExtSource), target :: external_source + ! Source and fission bank type(Neutron), allocatable, target :: source_bank(:) type(Bank), allocatable, target :: fission_bank(:) @@ -69,10 +70,10 @@ module global & SRC_SURFACE = 3 ! Source on a surface ! Problem type - integer :: problem_type = 0 + integer :: problem_type integer, parameter :: & - & PROB_CRITICALITY = 1, & ! Criticality problem - & PROB_SOURCE = 2 ! External source problem + & PROB_SOURCE = 1, & ! External source problem + & PROB_CRITICALITY = 2 ! Criticality problem ! Particle type integer, parameter :: & @@ -92,6 +93,20 @@ module global contains +!===================================================================== +! SET_DEFAULTS gives default values for many global parameters +!===================================================================== + + subroutine set_defaults() + + ! Default problem type is external source + problem_type = PROB_SOURCE + + ! Default number of particles + n_particles = 10000 + + end subroutine set_defaults + !===================================================================== ! FREE_MEMORY deallocates all allocatable arrays in the program, ! namely the cells, surfaces, materials, and sources diff --git a/src/input_sample b/src/input_sample index 80e83eda3..333ddb426 100644 --- a/src/input_sample +++ b/src/input_sample @@ -1,14 +1,14 @@ # test input file -cell 1 1 -2 1 -cell 2 1 -22 6 +cell 1 1 -2 +cell 2 1 2 -22 -surface 1 px 3.0 -surface 2 sph 0 0.0 0.0 5 +surface 1 sph 0 0 0 1 +surface 2 sph 0 0.0 0.0 3 material 1 U235 1.00 -source 1 box 0. 0. 0. 1. 1. 1. +source box -3 -3 -3 3 3 3 xs_library endfb7 -criticality 1000 100 1 +criticality 1 1 10 diff --git a/src/main.f90 b/src/main.f90 index 91071c760..c837d8438 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -5,7 +5,8 @@ program main use output, only: title, echo_input, message, warning, error use geometry, only: sense, cell_contains use mcnp_random, only: RN_init_problem, rang, RN_init_particle - use source, only: init_source + use source, only: init_source, get_source_particle + use physics, only: transport implicit none @@ -21,6 +22,9 @@ program main ! Initialize random number generator call RN_init_problem( 3, 0_8, 0_8, 0_8, 0 ) + ! Set default values for settings + call set_defaults() + ! Read command line arguments call read_command_line() @@ -35,7 +39,7 @@ program main call init_source() ! start problem - ! call run_problem() + call run_problem() contains @@ -59,14 +63,18 @@ contains HISTORY_LOOP: do j = 1, n_particles ! grab source particle from bank - ! particle => get_source_particle() + particle => get_source_particle() + if ( .not. associated(particle) ) then + ! no particles left in source bank + exit HISTORY_LOOP + end if ! set random number seed i_particle = (i-1)*n_particles + j call RN_init_particle(int(i_particle,8)) ! transport particle - ! call transport(particle) + call transport(particle) end do HISTORY_LOOP diff --git a/src/physics.f90 b/src/physics.f90 new file mode 100644 index 000000000..6130786df --- /dev/null +++ b/src/physics.f90 @@ -0,0 +1,57 @@ +module physics + + use global + use geometry, only: find_cell, dist_to_boundary + use types, only: Neutron + + implicit none + +contains + +!===================================================================== +! TRANSPORT encompasses the main logic for moving a particle through +! geometry. +!===================================================================== + + subroutine transport( neut ) + + type(Neutron), pointer, intent(inout) :: neut + + character(250) :: msg + integer :: i, surf + logical :: alive + + real(8) :: d_to_boundary + real(8) :: d_to_collision + real(8) :: distance + + if ( neut%cell == 0 ) then + call find_cell( neut ) + end if + + do while ( neut%alive ) + + ! Determine distance neutron moves + call dist_to_boundary( neut, d_to_boundary, surf ) + print *, d_to_boundary, surf + d_to_collision = -log(rang()) / 1.0 + distance = min( d_to_boundary, d_to_collision ) + + ! Advance neutron + neut%xyz = p%xyz + distance*p%uvw + + ! Add pathlength tallies + + if ( d_to_collision > d_to_boundary ) then + neut%surface = surf + neut%cell = 0 + call cross_boundary( neut ) + else + ! collision + end if + + end do + + end subroutine transport + +end module physics diff --git a/src/source.f90 b/src/source.f90 index 4d153fbe2..086e6e861 100644 --- a/src/source.f90 +++ b/src/source.f90 @@ -5,6 +5,8 @@ module source implicit none + integer :: source_index + contains !===================================================================== @@ -15,19 +17,68 @@ contains integer :: i,j real(8) :: r(3) + real(8) :: phi ! azimuthal angle + real(8) :: mu ! cosine of polar angle + real(8) :: p_min(3), p_max(3) ! Allocate fission and source banks allocate( source_bank(n_particles) ) allocate( fission_bank(3*n_particles) ) + ! Check external source type + if ( external_source%type == SRC_BOX ) then + p_min = external_source%values(1:3) + p_max = external_source%values(4:6) + end if + ! Initialize first cycle source bank do i = 1, n_particles call RN_init_particle(int(i,8)) + + ! position r = (/ (rang(), j = 1,3) /) source_bank(i)%uid = i - source_bank(i)%xyz = r + source_bank(i)%xyz = p_min + r*(p_max - p_min) + + ! angle + phi = 2.*PI*rang() + mu = 2.*rang() - 1 + source_bank(i)%uvw(1) = mu + source_bank(i)%uvw(2) = sqrt(1. - mu**2) * cos(phi) + source_bank(i)%uvw(3) = sqrt(1. - mu**2) * sin(phi) + + ! cell + source_bank(i)%cell = 0 + + ! set particle to be alive + source_bank(i)%alive = .true. end do + ! Reset source index + source_index = 0 + end subroutine init_source +!===================================================================== +! GET_SOURCE_PARTICLE returns the next source particle +!===================================================================== + + function get_source_particle() + + type(Neutron), pointer :: get_source_particle + + ! increment index + source_index = source_index + 1 + + ! if at end of bank, return null pointer + if ( source_index > size(source_bank) ) then + get_source_particle => null() + return + end if + + ! point to next source particle + get_source_particle => source_bank(source_index) + + end function get_source_particle + end module source diff --git a/src/types.f90 b/src/types.f90 index f553ad7d7..a66c55156 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -19,10 +19,13 @@ module types end type Cell type :: Neutron - integer :: uid ! Unique ID - real(8) :: xyz(3) ! location - real(8) :: uvw(3) ! directional cosines - real(8) :: wgt ! particle weight + integer :: uid ! Unique ID + real(8) :: xyz(3) ! location + real(8) :: uvw(3) ! directional cosines + integer :: cell ! current cell + integer :: surface ! current surface + real(8) :: wgt ! particle weight + logical :: alive ! is particle alive? end type Neutron type :: Bank @@ -44,7 +47,6 @@ module types end type Material type :: ExtSource - integer :: uid ! unique identifier integer :: type ! type of source, e.g. 'box' or 'cell' real(8), allocatable :: values(:) ! values for particular source type end type ExtSource