mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 21:25:36 -04:00
More constructive solid geometry routines.
This commit is contained in:
Paul Romano
parent
6d5cafbcb7
commit
55773c8e43
10 changed files with 347 additions and 108 deletions
13
ChangeLog
13
ChangeLog
|
|
@ -1,3 +1,16 @@
|
|||
2011-01-23 Paul Romano <romano7@mit.edu>
|
||||
|
||||
* 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 <romano7@mit.edu>
|
||||
|
||||
* input_sample: Added a sample input file
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
212
src/geometry.f90
212
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
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
16
src/main.f90
16
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
|
||||
|
||||
|
|
|
|||
57
src/physics.f90
Normal file
57
src/physics.f90
Normal file
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue