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

Initial commit with some basic geometry and IO routines.

Browse source
This commit is contained in:
Paul Romano 2011-01-17 23:52:55 +00:00
commit 03bbf502df
9 changed files with 1309 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

0
LICENSE Normal file
View file

0
README Normal file
View file

297
src/fileio.f90 Normal file
View file

@ -0,0 +1,297 @@
module fileio
use global
use types, only: Cell, Surface, CellList, SurfaceList
use string, only: split_string_wl, lower_case, string_to_real, split_string
use output, only: message, warning, error
implicit none
contains
!------------------------------------------------------------------------------
subroutine read_command_line()
integer :: argc
character(250) :: msg
logical :: file_exists
character(7) :: readable
argc = COMMAND_ARGUMENT_COUNT()
if ( argc > 0 ) then
call GET_COMMAND_ARGUMENT(1,inputfile)
else
msg = "No input file specified!"
call error( msg )
end if
! Check if input file exists
inquire( FILE=inputfile, EXIST=file_exists, READ=readable )
if ( .not. file_exists ) then
msg = "Input file '" // trim(inputfile) // "' does not exist!"
call error( msg )
elseif ( readable(1:3) /= 'YES' ) then
msg = "Input file '" // trim(inputfile) // "' is not readable! &
&Change file permissions with chmod command."
call error( msg )
end if
end subroutine read_command_line
!------------------------------------------------------------------------------
subroutine read_input(filename)
character(*), intent(in) :: filename
character(250) :: line
character(32) :: words(max_words)
integer :: in = 7
integer :: ioError
integer :: n, i
type(CellList), pointer :: c_head => null()
type(CellList), pointer :: c_tail => null()
type(CellList), pointer :: c_current => null()
type(SurfaceList), pointer :: s_head => null()
type(SurfaceList), pointer :: s_tail => null()
type(SurfaceList), pointer :: s_current => null()
!!$ type(MaterialList), pointer :: m_head
!!$ type(MaterialList), pointer :: m_tail
!!$ type(MaterialList), pointer :: m_current
ncell = 0
nsurf = 0
nmat = 0
open(file=filename, unit=in, status='old', action='read')
do
read(unit=in, fmt='(A250)', iostat=ioError) line
if ( ioError /= 0 ) exit
call split_string_wl(line, words, n)
select case ( trim(words(1)) )
case ('cell')
! Allocate a new cell
if ( .not. associated(c_head) ) then
allocate( c_head )
c_tail => c_head
else
allocate( c_tail%next )
c_tail => c_tail%next
end if
! Parse data
call parse_cell( c_tail, words, n )
ncell = ncell + 1
case ('surface')
! Allocate a new surface
if ( .not. associated(s_head) ) then
allocate( s_head )
s_tail => s_head
else
allocate( s_tail%next )
s_tail => s_tail%next
end if
! Parse data
call parse_surface( s_tail, words, n )
nsurf = nsurf + 1
case ( 'material' )
nmat = nmat + 1
case default
! do nothing
end select
end do
close(unit=in)
! Convert linked list of Cells into normal array
allocate( cells(ncell) )
c_current => c_head
do i = 1, ncell
! Allocate space on cells(i)
n = size(c_current%boundary_list)
allocate( cells(i)%boundary_list(n) )
! Copy cell
cells(i)%id = c_current%id
cells(i)%material = c_current%material
cells(i)%boundary_list = c_current%boundary_list
!!$ print *, 'Cell: ', cells(i)%id
!!$ print *, ' Material ', cells(i)%material
!!$ print *, ' Surfaces ', cells(i)%boundary_list
c_head => c_head%next
deallocate( c_current )
c_current => c_head
end do
! Convert linked list of Surfaces into normal array
allocate( surfaces(nsurf) )
s_current => s_head
do i = 1, nsurf
! Allocate space on surfaces(i)
n = size(s_current%coeffs)
allocate( surfaces(i)%coeffs(n) )
! Copy surface
surfaces(i)%id = s_current%id
surfaces(i)%type = s_current%type
surfaces(i)%coeffs = s_current%coeffs
!!$ print *, 'Surface: ', surfaces(i)%id
!!$ print *, ' Type ', surfaces(i)%type
!!$ print *, ' Coeffs ', surfaces(i)%coeffs
s_head => s_head%next
deallocate( s_current )
s_current => s_head
end do
! Convert linked list of Materials into normal array
end subroutine read_input
!------------------------------------------------------------------------------
subroutine parse_cell( cell, words, n_words )
type(CellList), pointer, intent(inout) :: cell
character(*), intent(in) :: words(max_words)
integer, intent(in) :: n_words
integer :: readError
integer :: i
character(250) :: msg
character(32) :: word
! Read cell identifier
read(words(2), fmt='(I8)', iostat=readError) cell%id
if ( readError > 0 ) then
msg = "Invalid cell name: " // words(2)
call error( msg )
end if
! Read cell material
read(words(3), fmt='(I8)', iostat=readError) cell%material
! Read list of surfaces
allocate( cell%boundary_list(n_words-3) )
do i = 1, n_words-3
word = words(i+3)
if ( word(1:1) == '(' ) then
cell%boundary_list(i) = OP_LEFT_PAREN
elseif ( word(1:1) == ')' ) then
cell%boundary_list(i) = OP_RIGHT_PAREN
elseif ( word(1:1) == ':' ) then
cell%boundary_list(i) = OP_UNION
elseif ( word(1:1) == '#' ) then
cell%boundary_list(i) = OP_DIFFERENCE
else
read(word, fmt='(I8)', iostat=readError) cell%boundary_list(i)
end if
end do
end subroutine parse_cell
!------------------------------------------------------------------------------
subroutine parse_surface( surface, words, n_words )
type(SurfaceList), pointer, intent(inout) :: surface
character(*), intent(in) :: words(max_words)
integer, intent(in) :: n_words
integer :: readError
integer :: i
integer :: coeffs_reqd
character(250) :: msg
character(32) :: word
! Read surface identifier
read(words(2), fmt='(I8)', iostat=readError) surface%id
if ( readError > 0 ) then
msg = "Invalid surface name: " // words(2)
call error( msg )
end if
! Read surface type
word = words(3)
call lower_case(word)
select case ( trim(word) )
case ( 'px' )
surface%type = SURF_PX
coeffs_reqd = 1
case ( 'py' )
surface%type = SURF_PY
coeffs_reqd = 1
case ( 'pz' )
surface%type = SURF_PZ
coeffs_reqd = 1
case ( 'plane' )
surface%type = SURF_PLANE
coeffs_reqd = 4
case ( 'cylx' )
surface%type = SURF_CYL_X
coeffs_reqd = 3
case ( 'cyly' )
surface%type = SURF_CYL_Y
coeffs_reqd = 3
case ( 'cylz' )
surface%type = SURF_CYL_Z
coeffs_reqd = 3
case ( 'sph' )
surface%type = SURF_SPHERE
coeffs_reqd = 4
case ( 'boxx' )
surface%type = SURF_BOX_X
coeffs_reqd = 4
case ( 'boxy' )
surface%type = SURF_BOX_Y
coeffs_reqd = 4
case ( 'boxz' )
surface%type = SURF_BOX_Z
coeffs_reqd = 4
case ( 'box' )
surface%type = SURF_BOX
coeffs_reqd = 6
case ( 'gq' )
surface%type = SURF_GQ
coeffs_reqd = 10
case default
msg = "Invalid surface type: " // words(3)
call error( msg )
end select
! Make sure there are enough coefficients for surface type
if ( n_words-3 < coeffs_reqd ) then
msg = "Not enough coefficients for surface: " // words(2)
call error( msg )
end if
! Read list of surfaces
allocate( surface%coeffs(n_words-3) )
do i = 1, n_words-3
surface%coeffs(i) = string_to_real(words(i+3))
end do
end subroutine parse_surface
!------------------------------------------------------------------------------
end module fileio

444
src/geometry.f90 Normal file
View file

@ -0,0 +1,444 @@
module geometry
use global
use types, only: Cell, Surface
use output, only: error
use string, only: int_to_string
implicit none
contains
!------------------------------------------------------------------------------
subroutine cell_contains( c, xyz, on_surface)
type(Cell), intent(in) :: c
real(8), intent(in) :: xyz(3)
logical, intent(in) :: on_surface
integer, allocatable :: expression(:)
integer :: specified_sense
integer :: actual_sense
integer :: n
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
! Lookup surface
surf_num = abs(expression(i))
! TODO: replace this loop with a hash since this lookup is O(N)
do j = 1,nsurf
surf => surfaces(j)
if ( surf%id == surf_num ) then
surface_found = .true.
exit
end if
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_string(surf_num)) // &
& " specified on cell " // int_to_string(c%id)
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
end if
! Need to deallocate expression
print *, expression
end subroutine cell_contains
!------------------------------------------------------------------------------
subroutine dist_to_boundary( cl, neut, dist )
type(Cell), intent(in) :: cl
type(Neutron), intent(in) :: neut
real(8), intent(out) :: dist
integer :: i
integer :: n_surfaces
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
real(8) :: tmp
real(8) :: a,b,c,k
real(8) :: quad
character(250) :: msg
x = neut%xyz(1)
y = neut%xyz(2)
z = neut%xyz(3)
u = neut%uvw(1)
v = neut%uvw(2)
z = neut%uvw(3)
dist = INFINITY
n_surfaces = size(cl%boundary_list)
do i = 1, n_surfaces
surf_num = abs(expression(i))
if ( surf_num >= OP_DIFFERENCE ) cycle
surf => surfaces(surf_num)
select case ( surf%type )
case ( SURF_PX )
if ( u == 0.0 ) then
d = INFINITY
else
x0 = surf%coeffs(1)
d = (x0 - x)/u
if ( d < 0 ) d = INFINITY
end if
case ( SURF_PY )
if ( v == 0.0 ) then
d = INFINITY
else
y0 = surf%coeffs(1)
d = (y0 - y)/v
if ( d < 0 ) d = INFINITY
end if
case ( SURF_PZ )
if ( w == 0.0 ) then
d = INFINITY
else
z0 = surf%coeffs(1)
d = (z0 - z)/w
if ( d < 0.0 ) d = INFINITY
end if
case ( SURF_PLANE )
A = surf%coeffs(1)
B = surf%coeffs(2)
C = surf%coeffs(3)
D = surf%coeffs(4)
tmp = A*u + B*v + C*w
if ( tmp == 0.0 ) then
d = INFINITY
else
d = -(A*x + B*y + C*w - D)/tmp
if ( d < 0.0 ) d = INFINITY
end if
case ( SURF_CYL_X )
a = 1.0 - u**2 ! v^2 + w^2
if ( a == 0.0 ) then
d = INFINITY
else
y0 = surf%coeffs(1)
z0 = surf%coeffs(2)
r = surf%coeffs(3)
y = y - y0
z = z - z0
k = y*v + z*w
c = y**2 + z**2 - r**2
quad = k**2 - a*c
if ( c < 0 ) then
! particle is inside the cylinder, thus one distance
! must be negative and one must be positive. The
! positive distance will be the one with negative sign
! on sqrt(quad)
d = -(k - sqrt(quad))/a
else
! particle is outside the cylinder, thus both
! distances are either positive or negative. If
! positive, the smaller distance is the one with
! positive sign on sqrt(quad)
d = -(k + sqrt(quad))/a
if ( d < 0 ) d = INFINITY
end if
end if
case ( SURF_CYL_Y )
a = 1.0 - v**2 ! u^2 + w^2
if ( a == 0.0 ) then
d = INFINITY
else
x0 = surf%coeffs(1)
z0 = surf%coeffs(2)
r = surf%coeffs(3)
x = x - x0
z = z - z0
k = x*u + z*w
c = x**2 + z**2 - r**2
quad = k**2 - a*c
if ( c < 0 ) then
! particle is inside the cylinder, thus one distance
! must be negative and one must be positive. The
! positive distance will be the one with negative sign
! on sqrt(quad)
d = -(k - sqrt(quad))/a
else
! particle is outside the cylinder, thus both
! distances are either positive or negative. If
! positive, the smaller distance is the one with
! positive sign on sqrt(quad)
d = -(k + sqrt(quad))/a
if ( d < 0 ) d = INFINITY
end if
end if
case ( SURF_CYL_Z )
a = 1.0 - w**2 ! u^2 + v^2
if ( a == 0.0 ) then
d = INFINITY
else
x0 = surf%coeffs(1)
y0 = surf%coeffs(2)
r = surf%coeffs(3)
x = x - x0
y = y - y0
k = x*u + y*v
c = x**2 + y**2 - r**2
quad = k**2 - a*c
if ( quad < 0 ) then
! no intersection with cylinder
d = INFINITY
elseif ( c < 0 ) then
! particle is inside the cylinder, thus one distance
! must be negative and one must be positive. The
! positive distance will be the one with negative sign
! on sqrt(quad)
d = -(k - sqrt(quad))/a
else
! particle is outside the cylinder, thus both
! distances are either positive or negative. If
! positive, the smaller distance is the one with
! positive sign on sqrt(quad)
d = -(k + sqrt(quad))/a
if ( d < 0 ) d = INFINITY
end if
end if
case ( SURF_SPHERE )
x0 = surf%coeffs(1)
y0 = surf%coeffs(2)
z0 = surf%coeffs(3)
r = surf%coeffs(4)
x = x - x0
y = y - y0
z = z - z0
k = x*u + y*v + z*w
c = x**2 + y**2 + z**2 - r**2
quad = k**2 - c
if ( quad < 0 ) then
! no intersection with sphere
d = INFINITY
elseif ( c < 0 ) then
! particle is inside the sphere, thus one distance
! must be negative and one must be positive. The
! positive distance will be the one with negative sign
! on sqrt(quad)
d = -(k - sqrt(quad))
else
! particle is outside the sphere, thus both
! distances are either positive or negative. If
! positive, the smaller distance is the one with
! positive sign on sqrt(quad)
d = -(k + sqrt(quad))
if ( d < 0 ) d = INFINITY
end if
case ( SURF_GQ )
msg = "Surface distance not yet implement for general quadratic."
call error( msg )
end select
! Check is calculated distance is new minimum
dist = min(d,dist)
end do
end subroutine dist_to_boundary
!------------------------------------------------------------------------------
subroutine sense( surf, xyz, s )
type(Surface), intent(in) :: surf
real(8), intent(in) :: xyz(3)
integer, intent(out) :: s
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
x = xyz(1)
y = xyz(2)
z = xyz(3)
select case ( surf%type )
case ( SURF_PX )
x0 = surf%coeffs(1)
func = x - x0
case ( SURF_PY )
y0 = surf%coeffs(1)
func = y - y0
case ( SURF_PZ )
z0 = surf%coeffs(1)
func = z - z0
case ( SURF_PLANE )
A = surf%coeffs(1)
B = surf%coeffs(2)
C = surf%coeffs(3)
D = surf%coeffs(4)
func = A*x + B*y + C*z - D
case ( SURF_CYL_X )
y0 = surf%coeffs(1)
z0 = surf%coeffs(2)
r = surf%coeffs(3)
func = (y-y0)**2 + (z-z0)**2 - r**2
case ( SURF_CYL_Y )
x0 = surf%coeffs(1)
z0 = surf%coeffs(2)
r = surf%coeffs(3)
func = (x-x0)**2 + (z-z0)**2 - r**2
case ( SURF_CYL_Z )
x0 = surf%coeffs(1)
y0 = surf%coeffs(2)
r = surf%coeffs(3)
func = (x-x0)**2 + (y-y0)**2 - r**2
case ( SURF_SPHERE )
x0 = surf%coeffs(1)
y0 = surf%coeffs(2)
z0 = surf%coeffs(3)
r = surf%coeffs(4)
func = (x-x0)**2 + (y-y0)**2 + (z-z0)**2 - r**2
case ( SURF_BOX_X )
y0 = surf%coeffs(1)
z0 = surf%coeffs(2)
y1 = surf%coeffs(3)
z1 = surf%coeffs(4)
if ( y >= y0 .and. y < y1 .and. z >= z0 .and. z < z1 ) then
s = SENSE_NEGATIVE
else
s = SENSE_POSITIVE
end if
return
case ( SURF_BOX_Y )
x0 = surf%coeffs(1)
z0 = surf%coeffs(2)
x1 = surf%coeffs(3)
z1 = surf%coeffs(4)
if ( x >= x0 .and. x < x1 .and. z >= z0 .and. z < z1 ) then
s = SENSE_NEGATIVE
else
s = SENSE_POSITIVE
end if
return
case ( SURF_BOX_Z )
x0 = surf%coeffs(1)
y0 = surf%coeffs(2)
x1 = surf%coeffs(3)
y1 = surf%coeffs(4)
if ( x >= x0 .and. x < x1 .and. y >= y0 .and. y < y1 ) then
s = SENSE_NEGATIVE
else
s = SENSE_POSITIVE
end if
return
case ( SURF_BOX )
x0 = surf%coeffs(1)
y0 = surf%coeffs(2)
z0 = surf%coeffs(3)
x1 = surf%coeffs(4)
y1 = surf%coeffs(5)
z1 = surf%coeffs(6)
if ( x >= x0 .and. x < x1 .and. y >= y0 .and. y < y1 .and. &
& z >= z0 .and. z < z1 ) then
s = SENSE_NEGATIVE
else
s = SENSE_POSITIVE
end if
return
case ( SURF_GQ )
func = A*x**2 + B*y**2 + C*z**2 + D*x*y + E*y*z + F*x*z + G*x &
& + H*y + I*z + J
end select
! Check which side of surface the point is on
if ( func > 0 ) then
s = SENSE_POSITIVE
else
s = SENSE_NEGATIVE
end if
end subroutine sense
end module geometry

88
src/global.f90 Normal file
View file

@ -0,0 +1,88 @@
module global
use types
implicit none
type(Cell), allocatable, target :: cells(:)
type(Surface), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
integer :: ncell ! Number of cells
integer :: nsurf ! Number of surfaces
integer :: nmat ! Number of materials
! Physical constants
real(8), parameter :: pi = 2.*acos(0.0) ! pi
! Boundary conditions
integer, parameter :: &
& TRANSMIT = 0, & ! Transmission boundary condition (default)
& VACUUM = 1, & ! Vacuum boundary condition
& REFLECT = 2 ! Reflecting boundary condition
! Logical operators for cell definitions
integer, parameter :: &
& OP_LEFT_PAREN = huge(0), & ! Left parentheses
& OP_RIGHT_PAREN = huge(0) - 1, & ! Right parentheses
& OP_UNION = huge(0) - 2, & ! Union operator
& OP_DIFFERENCE = huge(0) - 3 ! Difference operator
! Surface types
integer, parameter :: &
& SURF_PX = 1, & ! Plane parallel to x-plane
& SURF_PY = 2, & ! Plane parallel to y-plane
& SURF_PZ = 3, & ! Plane parallel to z-plane
& SURF_PLANE = 4, & ! Arbitrary plane
& SURF_CYL_X = 5, & ! Cylinder along x-axis
& SURF_CYL_Y = 6, & ! Cylinder along y-axis
& SURF_CYL_Z = 7, & ! Cylinder along z-axis
& SURF_SPHERE = 8, & ! Sphere
& SURF_BOX_X = 9, & ! Box extending infinity in x-direction
& SURF_BOX_Y = 10, & ! Box extending infinity in y-direction
& SURF_BOX_Z = 11, & ! Box extending infinity in z-direction
& SURF_BOX = 12, & ! Rectangular prism
& SURF_GQ = 13 ! General quadratic surface
integer, parameter :: &
& SENSE_POSITIVE = 1, &
& SENSE_NEGATIVE = -1
real(8), parameter :: &
& INFINITY = huge(0.0_8)
character(32) :: inputfile
integer, parameter :: verbosity = 5
integer, parameter :: max_words = 100
! Versioning numbers
integer, parameter :: VERSION_MAJOR = 0
integer, parameter :: VERSION_MINOR = 1
integer, parameter :: VERSION_RELEASE = 1
contains
!-----------------------------------------------------------------------
subroutine free_memory()
! Deallocate cells, surfaces, materials
if (allocated(cells)) then
deallocate(cells)
end if
if (allocated(surfaces)) then
deallocate(surfaces)
end if
if (allocated(materials)) then
deallocate(materials)
end if
! End program
stop
end subroutine free_memory
!-----------------------------------------------------------------------
end module global

26
src/main.f90 Normal file
View file

@ -0,0 +1,26 @@
program main
use global, only: cells, surfaces, materials, inputfile
use fileio, only: read_input, read_command_line
use output, only: title, message, warning, error
use geometry, only: sense, cell_contains
implicit none
character(16) :: filename
character(250) :: msg
real(8) :: point(3)
integer :: s(4)
call title()
call read_command_line()
call read_input(inputfile)
point = (/ 4.0, 0.0, 3.1 /)
call cell_contains( cells(1), point, .true. )
end program main

155
src/output.f90 Normal file
View file

@ -0,0 +1,155 @@
module output
use ISO_FORTRAN_ENV
use global, only: PATHLENGTH, COLLISION, ABSORPTION, FISSION, &
& NUFISSION, nhist, ncycle, localmesh, ijk_to_n, &
& nprocs, verbosity, VERSION_MAJOR, VERSION_MINOR, &
& VERSION_RELEASE, free_memory
implicit none
contains
!-----------------------------------------------------------------------
subroutine title()
character(10) :: date
character(8) :: time
write(OUTPUT_UNIT,*)
write(OUTPUT_UNIT,*) ' .d88888b. 888b d888 .d8888b. '
write(OUTPUT_UNIT,*) 'd88P" "Y88b 8888b d8888 d88P Y88b'
write(OUTPUT_UNIT,*) '888 888 88888b.d88888 888 888'
write(OUTPUT_UNIT,*) '888 888 88888b. .d88b. 88888b. 888Y88888P888 888 '
write(OUTPUT_UNIT,*) '888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 '
write(OUTPUT_UNIT,*) '888 888 888 888 88888888 888 888 888 Y8P 888 888 888'
write(OUTPUT_UNIT,*) 'Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P'
write(OUTPUT_UNIT,*) ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P" '
write(OUTPUT_UNIT,*) '____________888________________________________________________'
write(OUTPUT_UNIT,*) ' 888 '
write(OUTPUT_UNIT,*) ' 888 '
write(OUTPUT_UNIT,*)
! Write version information
! write(OUTPUT_UNIT,*) ' ______________________________________________________'
write(OUTPUT_UNIT,*) ' Developed At: Massachusetts Institute of Technology'
write(OUTPUT_UNIT,*) ' Lead Developer: Paul K. Romano'
write(OUTPUT_UNIT,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(OUTPUT_UNIT,101) trim(date), trim(time)
101 format (6X,"Date/Time:",7X,A,1X,A)
! write(OUTPUT_UNIT,*) ' ______________________________________________________'
write(OUTPUT_UNIT,*)
end subroutine title
!-----------------------------------------------------------------------
subroutine message( msg, level )
character(*), intent(in) :: msg
integer, intent(in) :: level
if ( level <= verbosity ) then
write (OUTPUT_UNIT,*) trim(msg)
end if
end subroutine message
!-----------------------------------------------------------------------
subroutine warning( msg )
character(*), intent(in) :: msg
integer :: n_lines
integer :: i
write(OUTPUT_UNIT, fmt='(1X,A9)', advance='no') 'WARNING: '
n_lines = (len_trim(msg)-1)/70 + 1
do i = 1, n_lines
if ( i == 1 ) then
write(OUTPUT_UNIT, fmt='(A70)') msg(70*(i-1)+1:70*i)
else
write(OUTPUT_UNIT, fmt='(10X,A70)') msg(70*(i-1)+1:70*i)
end if
end do
end subroutine warning
!-----------------------------------------------------------------------
subroutine error( msg )
character(*), intent(in) :: msg
integer :: n_lines
integer :: i
write(ERROR_UNIT, fmt='(1X,A7)', advance='no') 'ERROR: '
n_lines = (len_trim(msg)-1)/72 + 1
do i = 1, n_lines
if ( i == 1 ) then
write(ERROR_UNIT, fmt='(A72)') msg(72*(i-1)+1:72*i)
else
write(ERROR_UNIT, fmt='(7X,A72)') msg(72*(i-1)+1:72*i)
end if
end do
write(ERROR_UNIT,*)
call free_memory()
end subroutine error
!-----------------------------------------------------------------------
subroutine get_today( today_date, today_time )
character(10) :: today_date
character(8) :: today_time
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(7:7) // "/" // 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(7:7) // "/" // 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
!-----------------------------------------------------------------------
end module output

215
src/string.f90 Normal file
View file

@ -0,0 +1,215 @@
module string
use global, only: max_words
use output, only: error
implicit none
contains
!-------------------------------------------------------------------------------
subroutine split_string(string, words, n)
! SPLIT_STRING takes a sentence and splits it into separate words
! much the Python string.split() method.
!
! Arguments:
! string = input line
! words = array of words
! n = total number of words
character(*), intent(in) :: string
character(*), intent(out) :: words(max_words)
integer, intent(out) :: n
integer :: i ! current index
integer :: i_start ! starting index of word
integer :: i_end ! ending index of word
i_start = 0
i_end = 0
n = 0
do i = 1, len(string)
if ((i_start == 0) .and. (string(i:i) /= ' ')) then
i_start = i
end if
if (i_start > 0) then
if (string(i:i) == ' ') i_end = i - 1
if (i == len(string)) i_end = i
if (i_end > 0) then
n = n + 1
words(n) = string(i_start:i_end)
! reset indices
i_start = 0
i_end = 0
end if
end if
end do
end subroutine split_string
!-------------------------------------------------------------------------------
subroutine split_string_wl(string, words, n)
! SPLIT_STRING_WL takes a string that includes logical expressions
! for a list of bounding surfaces in a cell and splits it into
! separate words. The characters (, ), :, and # count as separate
! words since they represent operators.
!
! Arguments:
! string = input line
! words = array of words
! n = number of words
character(*), intent(in) :: string
character(*), intent(out) :: words(max_words)
integer, intent(out) :: n
character(1) :: char ! current character
integer :: i ! current index
integer :: i_start ! starting index of word
integer :: i_end ! ending index of word
i_start = 0
i_end = 0
n = 0
do i = 1, len_trim(string)
char = string(i:i)
! Check for special characters
if (index('():#', char) > 0) then
if (i_start > 0 ) then
i_end = i - 1
n = n + 1
words(n) = string(i_start:i_end)
end if
n = n + 1
words(n) = char
i_start = 0
i_end = 0
cycle
end if
if ((i_start == 0) .and. (char /= ' ')) then
i_start = i
end if
if (i_start > 0) then
if (char == ' ') i_end = i - 1
if (i == len_trim(string)) i_end = i
if (i_end > 0) then
n = n + 1
words(n) = string(i_start:i_end)
! reset indices
i_start = 0
i_end = 0
end if
end if
end do
end subroutine split_string_wl
!-------------------------------------------------------------------------------
subroutine concatenate( words, n_words, string )
! CONCATENATE takes an array of words and concatenates them
! together in one string with a single space between words
!
! Arguments:
! words = array of words
! n_words = total number of words
! string = concatenated string
character(*), intent(in) :: words(n_words)
integer, intent(in) :: n_words
character(250), intent(out) :: string
integer :: i ! index
string = words(1)
if ( n_words == 1 ) return
do i = 2, n_words
string = trim(string) // ' ' // words(i)
end do
end subroutine concatenate
!-------------------------------------------------------------------------------
elemental subroutine lower_case(word)
! convert a word to lower case
character(*), intent(inout) :: word
integer :: i
integer :: ic
do i = 1,len(word)
ic = ichar(word(i:i))
if (ic >= 65 .and. ic < 90) word(i:i) = char(ic+32)
end do
end subroutine lower_case
!-------------------------------------------------------------------------------
function string_to_real( string )
! STRING_TO_REAL converts an arbitrary string to a real(8)
!
! Arguments:
! string = character(*) containing number to convert
character(*), intent(in) :: string
real(8) :: string_to_real
integer :: index_decimal ! index of decimal point
integer :: index_exponent ! index of exponent character
integer :: w ! total field width
integer :: d ! number of digits to right of decimal point
integer :: readError
character(8) :: fmt ! format for reading string
character(250) :: msg ! error message
! Determine total field width
w = len_trim(string)
! Determine number of digits to right of decimal point
index_decimal = index(string, '.')
index_exponent = max(index(string, 'd'), index(string, 'D'), &
& index(string, 'e'), index(string, 'E'))
if ( index_decimal > 0 ) then
if ( index_exponent > 0 ) then
d = index_exponent - index_decimal - 1
else
d = w - index_decimal
end if
else
d = 0
end if
! Create format specifier for reading string
write( fmt, '("(E",I2,".",I2,")")') w, d
! Read string
read( string, fmt, iostat=readError ) string_to_real
if ( readError > 0 ) then
msg = "Could not read value: " // string
call error( msg )
end if
end function string_to_real
!-------------------------------------------------------------------------------
function int_to_string( num )
integer, intent(in) :: num
character(10) :: int_to_string
write ( int_to_string, '(I10)' ) num
int_to_string = adjustl(int_to_string)
end function int_to_string
!-------------------------------------------------------------------------------
end module string

84
src/types.f90 Normal file
View file

@ -0,0 +1,84 @@
module types
implicit none
type :: Material
real(8) :: Sigma_f ! Macro fission xs
real(8) :: Sigma_s ! Macro scattering xs
real(8) :: Sigma_a ! Macro absorption xs
real(8) :: Sigma_t ! Macro total xs
real(8) :: Sigma_nf ! Macro nu-fission xs
end type Material
type :: region
real(8) :: xyz_min(3) ! Minimum coordinates of region
real(8) :: xyz_max(3) ! Maximum coordinates of region
integer :: mat ! Material in region
end type region
type :: bank
real(8) :: xyz(3) ! location
integer :: uid ! Unique ID
end type bank
type :: particle
real(8) :: xyz(3) ! location
real(8) :: uvw(3) ! directional cosines
real(8) :: wgt ! particle weight
integer :: ijk(3) ! interval in global mesh
integer :: local_ijk(3) ! interval in local mesh
integer :: uid ! Unique ID
end type particle
type :: fissionBankPtr
type(bank), pointer :: ptr ! Used for sorting fission bank
end type fissionBankPtr
! For each basic type (Cell, Surface, and Material), we have
! corresponding types CellList, etc that have an extra pointer for a
! linked list. As the input file is being read, the list types are
! used since they can be dynamically allocated. Once the input file
! is complete, the linked list is converted to a normal array.
type :: Surface
integer :: id
integer :: type
real(8), allocatable :: coeffs(:)
integer, allocatable :: neighbor_pos(:)
integer, allocatable :: neighbor_neg(:)
integer :: bc
end type Surface
type :: SurfaceList
integer :: id
integer :: type
real(8), allocatable :: coeffs(:)
integer, allocatable :: neighbor_pos(:)
integer, allocatable :: neighbor_neg(:)
integer :: bc
type(SurfaceList), pointer :: next
end type SurfaceList
type :: Cell
integer :: id
integer :: n_items
integer, allocatable :: boundary_list(:)
integer :: material
end type Cell
type :: CellList
integer :: id
integer, allocatable :: boundary_list(:)
integer :: material
type(CellList), pointer :: next
end type CellList
type :: Neutron
integer :: uid ! Unique ID
real(8) :: xyz(3) ! location
real(8) :: uvw(3) ! directional cosines
real(8) :: wgt ! particle weight
end type Neutron
end module types