commit 03bbf502dfc4352d7e2c3cd7263a7c26ab050ee2 Author: Paul Romano Date: Mon Jan 17 23:52:55 2011 +0000 Initial commit with some basic geometry and IO routines. diff --git a/LICENSE b/LICENSE new file mode 100644 index 000000000..e69de29bb diff --git a/README b/README new file mode 100644 index 000000000..e69de29bb diff --git a/src/fileio.f90 b/src/fileio.f90 new file mode 100644 index 000000000..d1d5d8352 --- /dev/null +++ b/src/fileio.f90 @@ -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 + + + + diff --git a/src/geometry.f90 b/src/geometry.f90 new file mode 100644 index 000000000..a2b599c77 --- /dev/null +++ b/src/geometry.f90 @@ -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 + + diff --git a/src/global.f90 b/src/global.f90 new file mode 100644 index 000000000..bae0699da --- /dev/null +++ b/src/global.f90 @@ -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 diff --git a/src/main.f90 b/src/main.f90 new file mode 100644 index 000000000..b4d571dbf --- /dev/null +++ b/src/main.f90 @@ -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 diff --git a/src/output.f90 b/src/output.f90 new file mode 100644 index 000000000..31c0b2593 --- /dev/null +++ b/src/output.f90 @@ -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 + + + + diff --git a/src/string.f90 b/src/string.f90 new file mode 100644 index 000000000..83a4b44fd --- /dev/null +++ b/src/string.f90 @@ -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 diff --git a/src/types.f90 b/src/types.f90 new file mode 100644 index 000000000..b760cec1f --- /dev/null +++ b/src/types.f90 @@ -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