diff --git a/ChangeLog b/ChangeLog new file mode 100644 index 000000000..b4656841b --- /dev/null +++ b/ChangeLog @@ -0,0 +1,19 @@ +2011-01-22 Paul Romano + + * input_sample: Added a sample input file + * main.f90: Added run_problem subroutine + * mcnp_random.f90: Added random number generator based on MCNP's + linear congruential generator + * fileio.f90: Completely changed how cells, surfaces, and + materials are allocated. Now there is a first pass through the + input file made by read_count which allows the arrays to be + pre-allocated, then read_input actually parses the input data. + * geometry.f90: Only formatting changes + * output.f90: Added the echo_input subroutine, made formatting + changes + * string.f90: Moved int_to_str to the global module, changed + string splitting to recognize tab characters + * types.f90: Deleted List types, added Isotope, Material, and + Source + * source.f90: Added init_source routine which initializes source + particles diff --git a/src/Makefile b/src/Makefile index 9190f95ce..2e11fd2a0 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,7 +1,7 @@ program = openmc src = main.f90 types.f90 global.f90 fileio.f90 output.f90 \ - string.f90 geometry.f90 + string.f90 geometry.f90 mcnp_random.f90 source.f90 objects = $(src:.f90=.o) #-------------------------------------------------------------------- @@ -32,7 +32,9 @@ neat: types.o: global.o: types.o string.o: global.o output.o -output.o: global.o +output.o: global.o +source.o: global.o mcnp_random.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 +main.o: global.o fileio.o output.o geometry.o mcnp_random.o \ + source.o diff --git a/src/fileio.f90 b/src/fileio.f90 index d1d5d8352..0bb37f984 100644 --- a/src/fileio.f90 +++ b/src/fileio.f90 @@ -1,15 +1,17 @@ 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 types, only: Cell, Surface, ExtSource + use string, only: split_string_wl, lower_case, str_to_real, split_string use output, only: message, warning, error implicit none contains -!------------------------------------------------------------------------------ +!===================================================================== +! READ_COMMAND_LINE reads all parameters from the command line +!===================================================================== subroutine read_command_line() @@ -26,7 +28,7 @@ contains call error( msg ) end if - ! Check if input file exists + ! Check if input file exists and is readable inquire( FILE=inputfile, EXIST=file_exists, READ=readable ) if ( .not. file_exists ) then msg = "Input file '" // trim(inputfile) // "' does not exist!" @@ -39,7 +41,78 @@ contains end subroutine read_command_line -!------------------------------------------------------------------------------ +!===================================================================== +! READ_COUNT makes a first pass through the input file to determine +! how many cells, surfaces, materials, sources, etc there are in the +! problem. +!===================================================================== + + subroutine read_count(filename) + + character(*), intent(in) :: filename + + character(250) :: line, msg + character(32) :: words(max_words) + integer :: in = 7 + integer :: readError + integer :: n + + n_cells = 0 + n_surfaces = 0 + n_materials = 0 + n_sources = 0 + + open(file=filename, unit=in, status='old', action='read') + + do + read(unit=in, fmt='(A250)', iostat=readError) line + if ( readError /= 0 ) exit + call split_string(line, words, n) + if ( n == 0 ) cycle + + select case ( trim(words(1)) ) + case ( 'cell' ) + n_cells = n_cells + 1 + case ( 'surface' ) + 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 + if ( n_cells == 0 ) then + msg = "No cells specified!" + call error( msg ) + elseif ( n_surfaces == 0 ) then + msg = "No surfaces specified!" + call error( msg ) + 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) + + ! Allocate arrays for cells, surfaces, etc. + allocate( cells(n_cells) ) + allocate( surfaces(n_surfaces) ) + allocate( materials(n_materials) ) + allocate( sources(n_sources) ) + + end subroutine read_count + +!===================================================================== +! READ_INPUT takes a second pass through the input file and parses all +! the data on each line, calling the appropriate subroutine for each +! type of data item. +!===================================================================== subroutine read_input(filename) @@ -50,23 +123,15 @@ contains integer :: in = 7 integer :: ioError integer :: n, i + integer :: index_cell + integer :: index_surface + integer :: index_material + integer :: index_source - 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 - + index_cell = 0 + index_surface = 0 + index_material = 0 + index_source = 0 open(file=filename, unit=in, status='old', action='read') @@ -75,37 +140,29 @@ contains if ( ioError /= 0 ) exit call split_string_wl(line, words, n) + ! skip comments + if ( n==0 .or. words(1)(1:1) == '#' ) cycle + 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 + ! Read data from cell entry + index_cell = index_cell + 1 + call read_cell( index_cell, words, n ) 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 + ! Read data + index_surface = index_surface + 1 + call read_surface( index_surface, words, n ) case ( 'material' ) - nmat = nmat + 1 + + case ( 'source' ) + index_source = index_source + 1 + call read_source( index_source, words, n ) + + case ( 'criticality' ) + call read_criticality( words, n ) + case default ! do nothing end select @@ -114,102 +171,62 @@ contains 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 -!------------------------------------------------------------------------------ +!===================================================================== +! READ_CELL parses the data on a cell card. +!===================================================================== - subroutine parse_cell( cell, words, n_words ) + subroutine read_cell( index, words, n_words ) - type(CellList), pointer, intent(inout) :: cell + integer, intent(in) :: index character(*), intent(in) :: words(max_words) - integer, intent(in) :: n_words + integer, intent(in) :: n_words integer :: readError integer :: i character(250) :: msg character(32) :: word + type(Cell), pointer :: this_cell => null() + + this_cell => cells(index) ! Read cell identifier - read(words(2), fmt='(I8)', iostat=readError) cell%id + read(words(2), fmt='(I8)', iostat=readError) this_cell%uid 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(words(3), fmt='(I8)', iostat=readError) this_cell%material ! Read list of surfaces - allocate( cell%boundary_list(n_words-3) ) + allocate( this_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 + this_cell%boundary_list(i) = OP_LEFT_PAREN elseif ( word(1:1) == ')' ) then - cell%boundary_list(i) = OP_RIGHT_PAREN + this_cell%boundary_list(i) = OP_RIGHT_PAREN elseif ( word(1:1) == ':' ) then - cell%boundary_list(i) = OP_UNION + this_cell%boundary_list(i) = OP_UNION elseif ( word(1:1) == '#' ) then - cell%boundary_list(i) = OP_DIFFERENCE + this_cell%boundary_list(i) = OP_DIFFERENCE else - read(word, fmt='(I8)', iostat=readError) cell%boundary_list(i) + read(word, fmt='(I8)', iostat=readError) this_cell%boundary_list(i) end if end do - end subroutine parse_cell + end subroutine read_cell -!------------------------------------------------------------------------------ +!===================================================================== +! READ_SURFACE parses the data on a surface card. +!===================================================================== - subroutine parse_surface( surface, words, n_words ) + subroutine read_surface( index, words, n_words ) - type(SurfaceList), pointer, intent(inout) :: surface + integer, intent(in) :: index character(*), intent(in) :: words(max_words) integer, intent(in) :: n_words @@ -218,9 +235,12 @@ contains integer :: coeffs_reqd character(250) :: msg character(32) :: word + type(Surface), pointer :: this_surface => null() + + this_surface => surfaces(index) ! Read surface identifier - read(words(2), fmt='(I8)', iostat=readError) surface%id + read(words(2), fmt='(I8)', iostat=readError) this_surface%uid if ( readError > 0 ) then msg = "Invalid surface name: " // words(2) call error( msg ) @@ -231,43 +251,43 @@ contains call lower_case(word) select case ( trim(word) ) case ( 'px' ) - surface%type = SURF_PX + this_surface%type = SURF_PX coeffs_reqd = 1 case ( 'py' ) - surface%type = SURF_PY + this_surface%type = SURF_PY coeffs_reqd = 1 case ( 'pz' ) - surface%type = SURF_PZ + this_surface%type = SURF_PZ coeffs_reqd = 1 case ( 'plane' ) - surface%type = SURF_PLANE + this_surface%type = SURF_PLANE coeffs_reqd = 4 case ( 'cylx' ) - surface%type = SURF_CYL_X + this_surface%type = SURF_CYL_X coeffs_reqd = 3 case ( 'cyly' ) - surface%type = SURF_CYL_Y + this_surface%type = SURF_CYL_Y coeffs_reqd = 3 case ( 'cylz' ) - surface%type = SURF_CYL_Z + this_surface%type = SURF_CYL_Z coeffs_reqd = 3 case ( 'sph' ) - surface%type = SURF_SPHERE + this_surface%type = SURF_SPHERE coeffs_reqd = 4 case ( 'boxx' ) - surface%type = SURF_BOX_X + this_surface%type = SURF_BOX_X coeffs_reqd = 4 case ( 'boxy' ) - surface%type = SURF_BOX_Y + this_surface%type = SURF_BOX_Y coeffs_reqd = 4 case ( 'boxz' ) - surface%type = SURF_BOX_Z + this_surface%type = SURF_BOX_Z coeffs_reqd = 4 case ( 'box' ) - surface%type = SURF_BOX + this_surface%type = SURF_BOX coeffs_reqd = 6 case ( 'gq' ) - surface%type = SURF_GQ + this_surface%type = SURF_GQ coeffs_reqd = 10 case default msg = "Invalid surface type: " // words(3) @@ -281,14 +301,117 @@ contains end if ! Read list of surfaces - allocate( surface%coeffs(n_words-3) ) + allocate( this_surface%coeffs(n_words-3) ) do i = 1, n_words-3 - surface%coeffs(i) = string_to_real(words(i+3)) + this_surface%coeffs(i) = str_to_real(words(i+3)) end do - end subroutine parse_surface + end subroutine read_surface -!------------------------------------------------------------------------------ +!===================================================================== +! READ_SOURCE parses the data on a source card. +!===================================================================== + + subroutine read_source( index, 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) + call lower_case(word) + select case ( trim(word) ) + case ( 'box' ) + this_source%type = SRC_BOX + values_reqd = 6 + case default + msg = "Invalid source type: " // words(3) + 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) + 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)) + end do + + end subroutine read_source + +!===================================================================== +! READ_XS_LIBRARY parses the data on a xs_library card. This card +! specifies what cross section library should be used by default +!===================================================================== + + subroutine read_xs_library( words, n_words ) + + character(*), intent(in) :: words(max_words) + integer, intent(in) :: n_words + + end subroutine read_xs_library + +!===================================================================== +! READ_CRITICALITY parses the data on a criticality card. This card +! specifies that the problem at hand is a criticality calculation and +! gives the number of cycles, inactive cycles, and particles per +! cycle. +!===================================================================== + + subroutine read_criticality( words, n_words ) + + character(*), intent(in) :: words(max_words) + integer, intent(in) :: n_words + + character(250) :: msg + integer :: readError + + ! Set problem type to criticality + problem_type = PROB_CRITICALITY + + ! Read number of cycles + read(words(2), fmt='(I8)', iostat=readError) n_cycles + if ( readError > 0 ) then + msg = "Invalid number of cycles: " // words(2) + call error( msg ) + end if + + ! Read number of inactive cycles + read(words(3), fmt='(I8)', iostat=readError) n_inactive + if ( readError > 0 ) then + msg = "Invalid number of inactive cycles: " // words(2) + call error( msg ) + end if + + ! Read number of particles + read(words(4), fmt='(I8)', iostat=readError) n_particles + if ( readError > 0 ) then + msg = "Invalid number of particles: " // words(2) + call error( msg ) + end if + + end subroutine read_criticality end module fileio diff --git a/src/geometry.f90 b/src/geometry.f90 index a2b599c77..465a8982e 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -3,15 +3,16 @@ module geometry use global use types, only: Cell, Surface use output, only: error - use string, only: int_to_string implicit none contains -!------------------------------------------------------------------------------ +!===================================================================== +! CELL_CONTAINS determines whether a given point is inside a cell +!===================================================================== - subroutine cell_contains( c, xyz, on_surface) + subroutine cell_contains( c, xyz, on_surface ) type(Cell), intent(in) :: c real(8), intent(in) :: xyz(3) @@ -40,9 +41,9 @@ contains ! Lookup surface surf_num = abs(expression(i)) ! TODO: replace this loop with a hash since this lookup is O(N) - do j = 1,nsurf + do j = 1, n_surfaces surf => surfaces(j) - if ( surf%id == surf_num ) then + if ( surf%uid == surf_num ) then surface_found = .true. exit end if @@ -51,8 +52,8 @@ contains ! 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) + msg = "Count not find surface " // trim(int_to_str(surf_num)) // & + & " specified on cell " // int_to_str(c%uid) call error( msg ) end if @@ -73,7 +74,11 @@ contains end subroutine cell_contains -!------------------------------------------------------------------------------ +!===================================================================== +! DIST_TO_BOUNDARY calculates the distance to the nearest boundary of +! the cell 'cl' for a particle 'neut' traveling in a certain +! direction. +!===================================================================== subroutine dist_to_boundary( cl, neut, dist ) @@ -82,7 +87,7 @@ contains real(8), intent(out) :: dist integer :: i - integer :: n_surfaces + integer :: n_cell_surfaces integer, allocatable :: expression(:) integer :: surf_num type(Surface), pointer :: surf => null() @@ -102,8 +107,8 @@ contains z = neut%uvw(3) dist = INFINITY - n_surfaces = size(cl%boundary_list) - do i = 1, n_surfaces + n_cell_surfaces = size(cl%boundary_list) + do i = 1, n_cell_surfaces surf_num = abs(expression(i)) if ( surf_num >= OP_DIFFERENCE ) cycle @@ -310,7 +315,11 @@ contains end subroutine dist_to_boundary -!------------------------------------------------------------------------------ +!===================================================================== +! SENSE determines whether a point is on the 'positive' or 'negative' +! side of a surface. This routine is crucial for determining what cell +! a particular point is in. +!===================================================================== subroutine sense( surf, xyz, s ) diff --git a/src/global.f90 b/src/global.f90 index 2859a6616..9e59178c4 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -4,16 +4,29 @@ module global implicit none - type(Cell), allocatable, target :: cells(:) - type(Surface), allocatable, target :: surfaces(:) - type(Material), allocatable, target :: materials(:) + ! Main arrays for cells, surfaces, materials + 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 - integer :: ncell ! Number of cells - integer :: nsurf ! Number of surfaces - integer :: nmat ! Number of materials + ! 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 + + ! Source and fission bank + type(Neutron), allocatable, target :: source_bank(:) + type(Bank), allocatable, target :: fission_bank(:) ! Physical constants - real(8), parameter :: pi = 2.*acos(0.0) ! pi + real(8), parameter :: & + & PI = 2.*acos(0.0), & ! pi + & INFINITY = huge(0.0_8) ! positive infinity ! Boundary conditions integer, parameter :: & @@ -22,14 +35,14 @@ module global & REFLECT = 2 ! Reflecting boundary condition ! Logical operators for cell definitions - integer, parameter :: & + 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 :: & + 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 @@ -44,12 +57,28 @@ module global & SURF_BOX = 12, & ! Rectangular prism & SURF_GQ = 13 ! General quadratic surface - integer, parameter :: & + ! Surface senses + integer, parameter :: & & SENSE_POSITIVE = 1, & & SENSE_NEGATIVE = -1 - real(8), parameter :: & - & INFINITY = huge(0.0_8) + ! Source types + integer, parameter :: & + & SRC_BOX = 1, & ! Source in a rectangular prism + & SRC_CELL = 2, & ! Source in a cell + & SRC_SURFACE = 3 ! Source on a surface + + ! Problem type + integer :: problem_type = 0 + integer, parameter :: & + & PROB_CRITICALITY = 1, & ! Criticality problem + & PROB_SOURCE = 2 ! External source problem + + ! Particle type + integer, parameter :: & + & NEUTRON_ = 1, & + & PHOTON_ = 2, & + & ELECTRON_ = 3 character(32) :: inputfile @@ -63,7 +92,10 @@ module global contains -!----------------------------------------------------------------------- +!===================================================================== +! FREE_MEMORY deallocates all allocatable arrays in the program, +! namely the cells, surfaces, materials, and sources +!===================================================================== subroutine free_memory() @@ -78,11 +110,32 @@ contains deallocate(materials) end if + ! Deallocate fission and source bank + if (allocated(fission_bank)) then + deallocate(fission_bank) + end if + if (allocated(source_bank)) then + deallocate(source_bank) + end if + ! End program stop end subroutine free_memory -!----------------------------------------------------------------------- +!===================================================================== +! INT_TO_STR converts an integer to a string. Right now, it is limited +! to integers less than 10 billion. +!===================================================================== + + function int_to_str( num ) + + integer, intent(in) :: num + character(10) :: int_to_str + + write ( int_to_str, '(I10)' ) num + int_to_str = adjustl(int_to_str) + + end function int_to_str end module global diff --git a/src/input_sample b/src/input_sample new file mode 100644 index 000000000..80e83eda3 --- /dev/null +++ b/src/input_sample @@ -0,0 +1,14 @@ +# test input file + +cell 1 1 -2 1 +cell 2 1 -22 6 + +surface 1 px 3.0 +surface 2 sph 0 0.0 0.0 5 + +material 1 U235 1.00 + +source 1 box 0. 0. 0. 1. 1. 1. + +xs_library endfb7 +criticality 1000 100 1 diff --git a/src/main.f90 b/src/main.f90 index b4d571dbf..91071c760 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -1,9 +1,11 @@ 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 + use global + use fileio, only: read_input, read_command_line, read_count + 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 implicit none @@ -13,14 +15,75 @@ program main real(8) :: point(3) integer :: s(4) + ! Print the OpenMC title and version/date/time information call title() - + + ! Initialize random number generator + call RN_init_problem( 3, 0_8, 0_8, 0_8, 0 ) + + ! Read command line arguments call read_command_line() + + ! Read input file -- make a first pass through the file to count + ! cells, surfaces, etc in order to allocate arrays, then do a second + ! pass to actually read values + call read_count(inputfile) call read_input(inputfile) + call echo_input() - point = (/ 4.0, 0.0, 3.1 /) + ! create source particles + call init_source() - call cell_contains( cells(1), point, .true. ) + ! start problem + ! call run_problem() +contains + +!===================================================================== +! RUN_PROBLEM encompasses all the main logic where iterations are +! performed over the cycles and histories. +!===================================================================== + + subroutine run_problem() + + integer :: i, j + integer :: i_cycle + integer :: i_particle + type(Neutron), pointer :: particle => null() + + CYCLE_LOOP: do i_cycle = 1, n_cycles + + ! Set all tallies to zero + + HISTORY_LOOP: do j = 1, n_particles + + ! grab source particle from bank + ! particle => get_source_particle() + + ! set random number seed + i_particle = (i-1)*n_particles + j + call RN_init_particle(int(i_particle,8)) + + ! transport particle + ! call transport(particle) + + end do HISTORY_LOOP + + ! Collect results and statistics + + ! print cycle information + + end do CYCLE_LOOP + + ! Collect all tallies and print + + ! print run time + + call free_memory() + + end subroutine run_problem + end program main + + diff --git a/src/mcnp_random.f90 b/src/mcnp_random.f90 new file mode 100644 index 000000000..ba309cbdd --- /dev/null +++ b/src/mcnp_random.f90 @@ -0,0 +1,614 @@ + +module mcnp_random + !======================================================================= + ! Description: + ! mcnp_random.F90 -- random number generation routines + !======================================================================= + ! This module contains: + ! + ! * Constants for the RN generator, including initial RN seed for the + ! problem & the current RN seed + ! + ! * MCNP interface routines: + ! - random number function: rang() + ! - RN initialization for problem: RN_init_problem + ! - RN initialization for particle: RN_init_particle + ! - get info on RN parameters: RN_query + ! - get RN seed for n-th history: RN_query_first + ! - set new RN parameters: RN_set + ! - skip-ahead in the RN sequence: RN_skip_ahead + ! - Unit tests: RN_test_basic, RN_test_skip, RN_test_mixed + ! + ! * For interfacing with the rest of MCNP, arguments to/from these + ! routines will have types of I8 or I4. + ! Any args which are to hold random seeds, multipliers, + ! skip-distance will be type I8, so that 63 bits can be held without + ! truncation. + ! + ! Revisions: + ! * 10-04-2001 - F Brown, initial mcnp version + ! * 06-06-2002 - F Brown, mods for extended generators + ! * 12-21-2004 - F Brown, added 3 of LeCuyer's 63-bit mult. RNGs + ! * 01-29-2005 - J Sweezy, Modify to use mcnp modules prior to automatic + ! io unit numbers. + ! * 12-02-2005 - F Brown, mods for consistency with C version + !======================================================================= + + !------------------- + ! MCNP output units + !------------------- + integer, parameter :: iuo = 6 + integer, parameter :: jtty = 6 + + PRIVATE + !--------------------------------------------------- + ! Kinds for LONG INTEGERS (64-bit) & REAL*8 (64-bit) + !--------------------------------------------------- + integer, parameter :: R8 = selected_real_kind(15,307) + integer, parameter :: I8 = selected_int_kind(18) + + !----------------------------------- + ! Public functions and subroutines for this module + !----------------------------------- + PUBLIC :: rang + PUBLIC :: RN_init_problem + PUBLIC :: RN_init_particle + PUBLIC :: RN_set + PUBLIC :: RN_query + PUBLIC :: RN_query_first + PUBLIC :: RN_update_stats + PUBLIC :: RN_test_basic + PUBLIC :: RN_test_skip + PUBLIC :: RN_test_mixed + + !------------------------------------- + ! Constants for standard RN generators + !------------------------------------- + type :: RN_GEN + integer :: index + integer(I8) :: mult ! generator (multiplier) + integer(I8) :: add ! additive constant + integer :: log2mod ! log2 of modulus, must be <64 + integer(I8) :: stride ! stride for particle skip-ahead + integer(I8) :: initseed ! default seed for problem + character(len=8) :: name + end type RN_GEN + + ! parameters for standard generators + integer, parameter :: n_RN_GEN = 7 + type(RN_GEN), SAVE :: standard_generator(n_RN_GEN) + data standard_generator / & + & RN_GEN( 1, 19073486328125_I8, 0_I8, 48, 152917_I8, 19073486328125_I8 , 'mcnp std' ), & + & RN_GEN( 2, 9219741426499971445_I8, 1_I8, 63, 152917_I8, 1_I8, 'LEcuyer1' ), & + & RN_GEN( 3, 2806196910506780709_I8, 1_I8, 63, 152917_I8, 1_I8, 'LEcuyer2' ), & + & RN_GEN( 4, 3249286849523012805_I8, 1_I8, 63, 152917_I8, 1_I8, 'LEcuyer3' ), & + & RN_GEN( 5, 3512401965023503517_I8, 0_I8, 63, 152917_I8, 1_I8, 'LEcuyer4' ), & + & RN_GEN( 6, 2444805353187672469_I8, 0_I8, 63, 152917_I8, 1_I8, 'LEcuyer5' ), & + & RN_GEN( 7, 1987591058829310733_I8, 0_I8, 63, 152917_I8, 1_I8, 'LEcuyer6' ) & + & / + + !----------------------------------------------------------------- + ! * Linear multiplicative congruential RN algorithm: + ! + ! RN_SEED = RN_SEED*RN_MULT + RN_ADD mod RN_MOD + ! + ! * Default values listed below will be used, unless overridden + !----------------------------------------------------------------- + integer, SAVE :: RN_INDEX = 1 + integer(I8), SAVE :: RN_MULT = 19073486328125_I8 + integer(I8), SAVE :: RN_ADD = 0_I8 + integer, SAVE :: RN_BITS = 48 + integer(I8), SAVE :: RN_STRIDE = 152917_I8 + integer(I8), SAVE :: RN_SEED0 = 19073486328125_I8 + integer(I8), SAVE :: RN_MOD = 281474976710656_I8 + integer(I8), SAVE :: RN_MASK = 281474976710655_I8 + integer(I8), SAVE :: RN_PERIOD = 70368744177664_I8 + real(R8), SAVE :: RN_NORM = 1._R8 / 281474976710656._R8 + + !------------------------------------ + ! Private data for a single particle + !------------------------------------ + integer(I8) :: RN_SEED = 19073486328125_I8 ! current seed + integer(I8) :: RN_COUNT = 0_I8 ! current counter + integer(I8) :: RN_NPS = 0_I8 ! current particle number + + common /RN_THREAD/ RN_SEED, RN_COUNT, RN_NPS + save /RN_THREAD/ + !$OMP THREADprivate ( /RN_THREAD/ ) + + !------------------------------------------ + ! Shared data, to collect info on RN usage + !------------------------------------------ + integer(I8), SAVE :: RN_COUNT_TOTAL = 0 ! total RN count all particles + integer(I8), SAVE :: RN_COUNT_STRIDE = 0 ! count for stride exceeded + integer(I8), SAVE :: RN_COUNT_MAX = 0 ! max RN count all particles + integer(I8), SAVE :: RN_COUNT_MAX_NPS = 0 ! part index for max count + + !--------------------------------------------------------------------- + ! Reference data: Seeds for case of init.seed = 1, + ! Seed numbers for index 1-5, 123456-123460 + !--------------------------------------------------------------------- + integer(I8), dimension(10,n_RN_GEN) :: RN_CHECK + data RN_CHECK / & + ! ***** 1 ***** mcnp standard gen ***** + & 19073486328125_I8, 29763723208841_I8, 187205367447973_I8, & + & 131230026111313_I8, 264374031214925_I8, 260251000190209_I8, & + & 106001385730621_I8, 232883458246025_I8, 97934850615973_I8, & + & 163056893025873_I8, & + ! ***** 2 ***** + & 9219741426499971446_I8, 666764808255707375_I8, 4935109208453540924_I8, & + & 7076815037777023853_I8, 5594070487082964434_I8, 7069484152921594561_I8, & + & 8424485724631982902_I8, 19322398608391599_I8, 8639759691969673212_I8, & + & 8181315819375227437_I8, & + ! ***** 3 ***** + & 2806196910506780710_I8, 6924308458965941631_I8, 7093833571386932060_I8, & + & 4133560638274335821_I8, 678653069250352930_I8, 6431942287813238977_I8, & + & 4489310252323546086_I8, 2001863356968247359_I8, 966581798125502748_I8, & + & 1984113134431471885_I8, & + ! ***** 4 ***** + & 3249286849523012806_I8, 4366192626284999775_I8, 4334967208229239068_I8, & + & 6386614828577350285_I8, 6651454004113087106_I8, 2732760390316414145_I8, & + & 2067727651689204870_I8, 2707840203503213343_I8, 6009142246302485212_I8, & + & 6678916955629521741_I8, & + ! ***** 5 ***** + & 3512401965023503517_I8, 5461769869401032777_I8, 1468184805722937541_I8, & + & 5160872062372652241_I8, 6637647758174943277_I8, 794206257475890433_I8, & + & 4662153896835267997_I8, 6075201270501039433_I8, 889694366662031813_I8, & + & 7299299962545529297_I8, & + ! ***** 6 ***** + & 2444805353187672469_I8, 316616515307798713_I8, 4805819485453690029_I8, & + & 7073529708596135345_I8, 3727902566206144773_I8, 1142015043749161729_I8, & + & 8632479219692570773_I8, 2795453530630165433_I8, 5678973088636679085_I8, & + & 3491041423396061361_I8, & + ! ***** 7 ***** + & 1987591058829310733_I8, 5032889449041854121_I8, 4423612208294109589_I8, & + & 3020985922691845009_I8, 5159892747138367837_I8, 8387642107983542529_I8, & + & 8488178996095934477_I8, 708540881389133737_I8, 3643160883363532437_I8, & + & 4752976516470772881_I8 / + !--------------------------------------------------------------------- + +CONTAINS + + !------------------------------------------------------------------- + + function rang() + ! MCNP random number generator + ! + ! *************************************** + ! ***** modifies RN_SEED & RN_COUNT ***** + ! *************************************** + implicit none + real(R8) :: rang + + RN_SEED = iand( iand( RN_MULT*RN_SEED, RN_MASK) + RN_ADD, RN_MASK) + rang = RN_SEED * RN_NORM + RN_COUNT = RN_COUNT + 1 + + return + end function rang + + !------------------------------------------------------------------- + + function RN_skip_ahead( seed, skip ) + ! advance the seed "skip" RNs: seed*RN_MULT^n mod RN_MOD + implicit none + integer(I8) :: RN_skip_ahead + integer(I8), intent(in) :: seed, skip + integer(I8) :: nskip, gen, g, inc, c, gp, rn, seed_old + + seed_old = seed + ! add period till nskip>0 + nskip = skip + do while( nskip<0_I8 ) + if( RN_PERIOD>0_I8 ) then + nskip = nskip + RN_PERIOD + else + nskip = nskip + RN_MASK + nskip = nskip + 1_I8 + endif + enddo + + ! get gen=RN_MULT^n, in log2(n) ops, not n ops ! + nskip = iand( nskip, RN_MASK ) + gen = 1 + g = RN_MULT + inc = 0 + c = RN_ADD + do while( nskip>0_I8 ) + if( btest(nskip,0) ) then + gen = iand( gen*g, RN_MASK ) + inc = iand( inc*g, RN_MASK ) + inc = iand( inc+c, RN_MASK ) + endif + gp = iand( g+1, RN_MASK ) + g = iand( g*g, RN_MASK ) + c = iand( gp*c, RN_MASK ) + nskip = ishft( nskip, -1 ) + enddo + rn = iand( gen*seed_old, RN_MASK ) + rn = iand( rn + inc, RN_MASK ) + RN_skip_ahead = rn + return + end function RN_skip_ahead + + !------------------------------------------------------------------- + + subroutine RN_init_problem( new_standard_gen, new_seed, & + & new_stride, new_part1, print_info ) + ! * initialize MCNP random number parameters for problem, + ! based on user input. This routine should be called + ! only from the main thread, if OMP threading is being used. + ! + ! * for initial & continue runs, these args should be set: + ! new_standard_gen - index of built-in standard RN generator, + ! from RAND gen= (or dbcn(14) + ! new_seed - from RAND seed= (or dbcn(1)) + ! output - logical, print RN seed & mult if true + ! + ! new_stride - from RAND stride= (or dbcn(13)) + ! new_part1 - from RAND hist= (or dbcn(8)) + ! + ! * for continue runs only, these should also be set: + ! new_count_total - from "rnr" at end of previous run + ! new_count_stride - from nrnh(1) at end of previous run + ! new_count_max - from nrnh(2) at end of previous run + ! new_count_max_nps - from nrnh(3) at end of previous run + ! + ! * check on size of long-ints & long-int arithmetic + ! * check the multiplier + ! * advance the base seed for the problem + ! * set the initial particle seed + ! * initialize the counters for RN stats + implicit none + integer, intent(in) :: new_standard_gen + integer(I8), intent(in) :: new_seed + integer(I8), intent(in) :: new_stride + integer(I8), intent(in) :: new_part1 + integer, intent(in) :: print_info + character(len=20) :: printseed + integer(I8) :: itemp1, itemp2, itemp3, itemp4 + + if( new_standard_gen<1 .or. new_standard_gen>n_RN_GEN ) then + call expire( 0, 'RN_init_problem', & + & ' ***** ERROR: illegal index for built-in RN generator') + endif + + ! set defaults, override if input supplied: seed, mult, stride + RN_INDEX = new_standard_gen + RN_MULT = standard_generator(RN_INDEX)%mult + RN_ADD = standard_generator(RN_INDEX)%add + RN_STRIDE = standard_generator(RN_INDEX)%stride + RN_SEED0 = standard_generator(RN_INDEX)%initseed + RN_BITS = standard_generator(RN_INDEX)%log2mod + RN_MOD = ishft( 1_I8, RN_BITS ) + RN_MASK = ishft( not(0_I8), RN_BITS-64 ) + RN_NORM = 2._R8**(-RN_BITS) + if( RN_ADD==0_I8) then + RN_PERIOD = ishft( 1_I8, RN_BITS-2 ) + else + RN_PERIOD = ishft( 1_I8, RN_BITS ) + endif + if( new_seed>0_I8 ) then + RN_SEED0 = new_seed + endif + if( new_stride>0_I8 ) then + RN_STRIDE = new_stride + endif + RN_COUNT_TOTAL = 0 + RN_COUNT_STRIDE = 0 + RN_COUNT_MAX = 0 + RN_COUNT_MAX_NPS = 0 + + if( print_info /= 0 ) then + write(printseed,'(i20)') RN_SEED0 + write( iuo,1) RN_INDEX, RN_SEED0, RN_MULT, RN_ADD, RN_BITS, RN_STRIDE + write(jtty,2) RN_INDEX, adjustl(printseed) +1 format( & + & /,' ***************************************************', & + & /,' * Random Number Generator = ',i20, ' *', & + & /,' * Random Number Seed = ',i20, ' *', & + & /,' * Random Number Multiplier = ',i20, ' *', & + & /,' * Random Number Adder = ',i20, ' *', & + & /,' * Random Number Bits Used = ',i20, ' *', & + & /,' * Random Number Stride = ',i20, ' *', & + & /,' ***************************************************',/) +2 format(' comment. using random number generator ',i2, & + & ', initial seed = ',a20) + endif + + ! double-check on number of bits in a long int + if( bit_size(RN_SEED)<64 ) then + call expire( 0, 'RN_init_problem', & + & ' ***** ERROR: <64 bits in long-int, can-t generate RN-s') + endif + itemp1 = 5_I8**25 + itemp2 = 5_I8**19 + itemp3 = ishft(2_I8**62-1_I8,1) + 1_I8 + itemp4 = itemp1*itemp2 + if( iand(itemp4,itemp3)/=8443747864978395601_I8 ) then + call expire( 0, 'RN_init_problem', & + & ' ***** ERROR: can-t do 64-bit integer ops for RN-s') + endif + + if( new_part1>1_I8 ) then + ! advance the problem seed to that for part1 + RN_SEED0 = RN_skip_ahead( RN_SEED0, (new_part1-1_I8)*RN_STRIDE ) + itemp1 = RN_skip_ahead( RN_SEED0, RN_STRIDE ) + if( print_info /= 0 ) then + write(printseed,'(i20)') itemp1 + write( iuo,3) new_part1, RN_SEED0, itemp1 + write(jtty,4) new_part1, adjustl(printseed) +3 format( & + & /,' ***************************************************', & + & /,' * Random Number Seed will be advanced to that for *', & + & /,' * previous particle number = ',i20, ' *', & + & /,' * New RN Seed for problem = ',i20, ' *', & + & /,' * Next Random Number Seed = ',i20, ' *', & + & /,' ***************************************************',/) +4 format(' comment. advancing random number to particle ',i12, & + & ', initial seed = ',a20) + endif + endif + + ! set the initial particle seed + RN_SEED = RN_SEED0 + RN_COUNT = 0 + RN_NPS = 0 + + return + end subroutine RN_init_problem + + !------------------------------------------------------------------- + + subroutine RN_init_particle( nps ) + ! initialize MCNP random number parameters for particle "nps" + ! + ! * generate a new particle seed from the base seed + ! & particle index + ! * set the RN count to zero + implicit none + integer(I8), intent(in) :: nps + + RN_SEED = RN_skip_ahead( RN_SEED0, nps*RN_STRIDE ) + RN_COUNT = 0 + RN_NPS = nps + + return + end subroutine RN_init_particle + + !------------------------------------------------------------------- + + subroutine RN_set( key, value ) + implicit none + character(len=*), intent(in) :: key + integer(I8), intent(in) :: value + character(len=20) :: printseed + integer(I8) :: itemp1 + + if( key == "stride" ) then + if( value>0_I8 ) then + RN_STRIDE = value + endif + endif + if( key == "count_total" ) RN_COUNT_TOTAL = value + if( key == "count_stride" ) RN_COUNT_STRIDE = value + if( key == "count_max" ) RN_COUNT_MAX = value + if( key == "count_max_nps" ) RN_COUNT_MAX_NPS = value + if( key == "seed" ) then + if( value>0_I8 ) then + RN_SEED0 = value + RN_SEED = RN_SEED0 + RN_COUNT = 0 + RN_NPS = 0 + endif + endif + if( key == "part1" ) then + if( value>1_I8 ) then + ! advance the problem seed to that for part1 + RN_SEED0 = RN_skip_ahead( RN_SEED0, (value-1_I8)*RN_STRIDE ) + itemp1 = RN_skip_ahead( RN_SEED0, RN_STRIDE ) + write(printseed,'(i20)') itemp1 + write( iuo,3) value, RN_SEED0, itemp1 + write(jtty,4) value, adjustl(printseed) +3 format( & + & /,' ***************************************************', & + & /,' * Random Number Seed will be advanced to that for *', & + & /,' * previous particle number = ',i20, ' *', & + & /,' * New RN Seed for problem = ',i20, ' *', & + & /,' * Next Random Number Seed = ',i20, ' *', & + & /,' ***************************************************',/) +4 format(' comment. advancing random number to particle ',i12, & + & ', initial seed = ',a20) + RN_SEED = RN_SEED0 + RN_COUNT = 0 + RN_NPS = 0 + endif + endif + return + end subroutine RN_set + + !------------------------------------------------------------------- + + function RN_query( key ) + implicit none + integer(I8) :: RN_query + character(len=*), intent(in) :: key + RN_query = 0_I8 + if( key == "seed" ) RN_query = RN_SEED + if( key == "stride" ) RN_query = RN_STRIDE + if( key == "mult" ) RN_query = RN_MULT + if( key == "add" ) RN_query = RN_ADD + if( key == "count" ) RN_query = RN_COUNT + if( key == "period" ) RN_query = RN_PERIOD + if( key == "count_total" ) RN_query = RN_COUNT_TOTAL + if( key == "count_stride" ) RN_query = RN_COUNT_STRIDE + if( key == "count_max" ) RN_query = RN_COUNT_MAX + if( key == "count_max_nps" ) RN_query = RN_COUNT_MAX_NPS + if( key == "first" ) RN_query = RN_SEED0 + return + end function RN_query + !------------------------------------------------------------------- + + function RN_query_first( nps ) + implicit none + integer(I8) :: RN_query_first + integer(I8), intent(in) :: nps + RN_query_first = RN_skip_ahead( RN_SEED0, nps*RN_STRIDE ) + return + end function RN_query_first + + !------------------------------------------------------------------- + + subroutine RN_update_stats() + ! update overall RN count info + implicit none + + !$OMP CRITICAL (RN_STATS) + + RN_COUNT_TOTAL = RN_COUNT_TOTAL + RN_COUNT + + if( RN_COUNT>RN_COUNT_MAX ) then + RN_COUNT_MAX = RN_COUNT + RN_COUNT_MAX_NPS = RN_NPS + endif + + if( RN_COUNT>RN_STRIDE ) then + RN_COUNT_STRIDE = RN_COUNT_STRIDE + 1 + endif + + !$OMP END CRITICAL (RN_STATS) + + RN_COUNT = 0 + RN_NPS = 0 + + return + end subroutine RN_update_stats + + !------------------------------------------------------------------- + + subroutine expire( i, c1, c2 ) + integer, intent(in) :: i + character(len=*), intent(in) :: c1, c2 + write(*,*) ' ********** error: ',c1 + write(*,*) ' ********** error: ',c2 + stop '**error**' + end subroutine expire + + !------------------------------------------------------------------- + !################################################################### + !# + !# Unit tests + !# + !################################################################### + + subroutine RN_test_basic( new_gen ) + ! test routine for basic random number generator + implicit none + integer, intent(in) :: new_gen + real(R8) :: s + integer(I8) :: seeds(10) + integer :: i, j + + write(jtty,"(/,a)") " ***** random number - basic test *****" + + ! set the seed + call RN_init_problem( new_gen, 1_I8, 0_I8, 0_I8, 1 ) + + ! get the first 5 seeds, then skip a few, get 5 more - directly + s = 0.0_R8 + do i = 1,5 + s = s + rang() + seeds(i) = RN_query( "seed" ) + enddo + do i = 6,123455 + s = s + rang() + enddo + do i = 6,10 + s = s + rang() + seeds(i) = RN_query( "seed" ) + enddo + + ! compare + do i = 1,10 + j = i + if( i>5 ) j = i + 123450 + write(jtty,"(1x,i6,a,i20,a,i20)") & + & j, " reference: ", RN_CHECK(i,new_gen), " computed: ", seeds(i) + if( seeds(i)/=RN_CHECK(i,new_gen) ) then + write(jtty,"(a)") " ***** basic_test of RN generator failed:" + endif + enddo + return + end subroutine RN_test_basic + + !------------------------------------------------------------------- + + subroutine RN_test_skip( new_gen ) + ! test routine for basic random number generation & skip-ahead + implicit none + integer, intent(in) :: new_gen + integer(I8) :: seeds(10) + integer :: i, j + + ! set the seed + call RN_init_problem( new_gen, 1_I8, 0_I8, 0_I8, 0 ) + + ! use the skip-ahead function to get first 5 seeds, then 5 more + do i = 1,10 + j = i + if( i>5 ) j = i + 123450 + seeds(i) = RN_skip_ahead( 1_I8, int(j,I8) ) + enddo + + ! compare + write(jtty,"(/,a)") " ***** random number - skip test *****" + do i = 1,10 + j = i + if( i>5 ) j = i + 123450 + write(jtty,"(1x,i6,a,i20,a,i20)") & + & j, " reference: ", RN_CHECK(i,new_gen), " computed: ", seeds(i) + if( seeds(i)/=RN_CHECK(i,new_gen) ) then + write(jtty,"(a)") " ***** skip_test of RN generator failed:" + endif + enddo + return + end subroutine RN_test_skip + + !------------------------------------------------------------------- + + subroutine RN_test_mixed( new_gen ) + ! test routine -- print RN's 1-5 & 123456-123460, + ! with reference vals + implicit none + integer, intent(in) :: new_gen + integer(I8) :: r + integer :: i, j + + write(jtty,"(/,a)") " ***** random number - mixed test *****" + ! set the seed & set the stride to 1 + call RN_init_problem( new_gen, 1_I8, 1_I8, 0_I8, 0 ) + + write(jtty,"(a,i20,z20)") " RN_MULT = ", RN_MULT, RN_MULT + write(jtty,"(a,i20,z20)") " RN_ADD = ", RN_ADD, RN_ADD + write(jtty,"(a,i20,z20)") " RN_MOD = ", RN_MOD, RN_MOD + write(jtty,"(a,i20,z20)") " RN_MASK = ", RN_MASK, RN_MASK + write(jtty,"(a,i20)") " RN_BITS = ", RN_BITS + write(jtty,"(a,i20)") " RN_PERIOD = ", RN_PERIOD + write(jtty,"(a,es20.14)") " RN_NORM = ", RN_NORM + write(jtty,"(a)") " " + do i = 1,10 + j = i + if( i>5 ) j = i + 123450 + call RN_init_particle( int(j,I8) ) + r = RN_query( "seed" ) + write(jtty,"(1x,i6,a,i20,a,i20)") & + & j, " reference: ", RN_CHECK(i,new_gen)," computed: ", r + if( r/=RN_CHECK(i,new_gen) ) then + write(jtty,"(a)") " ***** mixed test of RN generator failed:" + endif + enddo + return + end subroutine RN_test_mixed + + !------------------------------------------------------------------- +end module mcnp_random diff --git a/src/output.f90 b/src/output.f90 index 425e70b15..d84da2697 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -1,51 +1,101 @@ module output use ISO_FORTRAN_ENV - use global, only: verbosity, VERSION_MAJOR, VERSION_MINOR, & - & VERSION_RELEASE, free_memory + use global implicit none contains -!----------------------------------------------------------------------- +!===================================================================== +! TITLE prints the main title banner as well as information about the +! program developers, version, and date/time which the problem was +! run. +!===================================================================== subroutine title() character(10) :: date character(8) :: time + integer :: ou - 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,*) + ou = OUTPUT_UNIT + + write(ou,*) + write(ou,*) ' .d88888b. 888b d888 .d8888b. ' + write(ou,*) 'd88P" "Y88b 8888b d8888 d88P Y88b' + write(ou,*) '888 888 88888b.d88888 888 888' + write(ou,*) '888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ' + write(ou,*) '888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ' + write(ou,*) '888 888 888 888 88888888 888 888 888 Y8P 888 888 888' + write(ou,*) 'Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P' + write(ou,*) ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P" ' + write(ou,*) '____________888________________________________________________' + write(ou,*) ' 888 ' + write(ou,*) ' 888 ' + write(ou,*) ! 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 +! write(ou,*) ' ______________________________________________________' + write(ou,*) ' Developed At: Massachusetts Institute of Technology' + write(ou,*) ' Lead Developer: Paul K. Romano' + write(ou,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) + write(ou,101) trim(date), trim(time) 101 format (6X,"Date/Time:",7X,A,1X,A) -! write(OUTPUT_UNIT,*) ' ______________________________________________________' - write(OUTPUT_UNIT,*) +! write(ou,*) ' ______________________________________________________' + write(ou,*) end subroutine title -!----------------------------------------------------------------------- +!===================================================================== +! ECHO_INPUT displays summary information about the problem about to +! be run after reading all the input. +!===================================================================== + + subroutine echo_input() + + character(32) :: string + integer :: ou + + ou = OUTPUT_UNIT + + ! Display problem summary + write(ou,*) '=============================================' + write(ou,*) '=> PROBLEM SUMMARY <=' + write(ou,*) '=============================================' + write(ou,*) + if ( problem_type == PROB_CRITICALITY ) then + write(ou,100) 'Problem type:', 'Criticality' + write(ou,100) 'Number of Cycles:', int_to_str(n_cycles) + write(ou,100) 'Number of Inactive Cycles:', int_to_str(n_inactive) + elseif ( problem_type == PROB_SOURCE ) then + write(ou,100) 'Problem type:', 'External Source' + end if + write(ou,100) 'Number of Particles:', int_to_str(n_particles) + write(ou,*) + ! Display geometry summary + write(ou,*) '=============================================' + write(ou,*) '=> GEOMETRY SUMMARY <=' + write(ou,*) '=============================================' + write(ou,*) + write(ou,100) 'Number of Cells:', int_to_str(n_cells) + write(ou,100) 'Number of Surfaces:', int_to_str(n_surfaces) + write(ou,100) 'Number of Materials:', int_to_str(n_materials) + write(ou,*) + + ! Format descriptor for columns +100 format (1X,A,T35,A) + + end subroutine echo_input + +!===================================================================== +! MESSAGE displays an informational message to the log file and the +! standard output stream. +!===================================================================== subroutine message( msg, level ) @@ -58,59 +108,73 @@ module output end subroutine message -!----------------------------------------------------------------------- +!===================================================================== +! WARNING issues a warning to the user in the log file and the +! standard output stream. +!===================================================================== subroutine warning( msg ) character(*), intent(in) :: msg integer :: n_lines - integer :: i + integer :: i, ou - write(OUTPUT_UNIT, fmt='(1X,A9)', advance='no') 'WARNING: ' + ou = OUTPUT_UNIT + + write(ou, 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) + write(ou, fmt='(A70)') msg(70*(i-1)+1:70*i) else - write(OUTPUT_UNIT, fmt='(10X,A70)') msg(70*(i-1)+1:70*i) + write(ou, fmt='(10X,A70)') msg(70*(i-1)+1:70*i) end if end do end subroutine warning -!----------------------------------------------------------------------- +!===================================================================== +! ERROR alerts the user that an error has been encountered and +! displays a message about the particular problem. Errors are +! considered 'fatal' and hence the program is aborted. +!===================================================================== subroutine error( msg ) character(*), intent(in) :: msg integer :: n_lines - integer :: i + integer :: i, eu - write(ERROR_UNIT, fmt='(1X,A7)', advance='no') 'ERROR: ' + eu = ERROR_UNIT + + write(eu, 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) + write(eu, fmt='(A72)') msg(72*(i-1)+1:72*i) else - write(ERROR_UNIT, fmt='(7X,A72)') msg(72*(i-1)+1:72*i) + write(eu, fmt='(7X,A72)') msg(72*(i-1)+1:72*i) end if end do - write(ERROR_UNIT,*) + write(eu,*) call free_memory() end subroutine error -!----------------------------------------------------------------------- +!===================================================================== +! GET_TODAY determines the date and time at which the program began +! execution and returns it in a readable format +!===================================================================== subroutine get_today( today_date, today_time ) - character(10) :: today_date - character(8) :: today_time + character(10), intent(out) :: today_date + character(8), intent(out) :: today_time character(8) :: date character(10) :: time @@ -144,10 +208,4 @@ module output end subroutine get_today -!----------------------------------------------------------------------- - end module output - - - - diff --git a/src/source.f90 b/src/source.f90 new file mode 100644 index 000000000..4d153fbe2 --- /dev/null +++ b/src/source.f90 @@ -0,0 +1,33 @@ +module source + + use global + use mcnp_random, only: rang, RN_init_particle + + implicit none + +contains + +!===================================================================== +! INIT_SOURCE initializes particles in the source bank +!===================================================================== + + subroutine init_source() + + integer :: i,j + real(8) :: r(3) + + ! Allocate fission and source banks + allocate( source_bank(n_particles) ) + allocate( fission_bank(3*n_particles) ) + + ! Initialize first cycle source bank + do i = 1, n_particles + call RN_init_particle(int(i,8)) + r = (/ (rang(), j = 1,3) /) + source_bank(i)%uid = i + source_bank(i)%xyz = r + end do + + end subroutine init_source + +end module source diff --git a/src/string.f90 b/src/string.f90 index 83a4b44fd..d0f391123 100644 --- a/src/string.f90 +++ b/src/string.f90 @@ -7,21 +7,23 @@ module string contains -!------------------------------------------------------------------------------- +!===================================================================== +! 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 +!===================================================================== 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 + character(1) :: char ! current character integer :: i ! current index integer :: i_start ! starting index of word integer :: i_end ! ending index of word @@ -30,11 +32,14 @@ contains i_end = 0 n = 0 do i = 1, len(string) - if ((i_start == 0) .and. (string(i:i) /= ' ')) then + char = string(i:i) + + ! Note that ACHAR(9) is a horizontal tab + if ((i_start == 0) .and. (char /= ' ') .and. (char /= achar(9))) then i_start = i end if if (i_start > 0) then - if (string(i:i) == ' ') i_end = i - 1 + if ((char == ' ') .or. (char == achar(9))) i_end = i - 1 if (i == len(string)) i_end = i if (i_end > 0) then n = n + 1 @@ -48,18 +53,19 @@ contains end subroutine split_string -!------------------------------------------------------------------------------- +!===================================================================== +! 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 +!===================================================================== 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) @@ -107,16 +113,17 @@ contains end do end subroutine split_string_wl -!------------------------------------------------------------------------------- +!===================================================================== +! 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 +!===================================================================== 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 @@ -132,7 +139,9 @@ contains end subroutine concatenate -!------------------------------------------------------------------------------- +!===================================================================== +! LOWER_CASE converts a string to all lower case characters +!===================================================================== elemental subroutine lower_case(word) ! convert a word to lower case @@ -149,16 +158,21 @@ contains end subroutine lower_case -!------------------------------------------------------------------------------- +!===================================================================== +! STR_TO_REAL converts an arbitrary string to a real(8). Generally +! this function is intended for strings for which the exact format is +! not known. If the format of the number is known a priori, the +! appropriate format descriptor should be used in lieu of this routine +! because of the extra overhead. +! +! Arguments: +! string = character(*) containing number to convert +!===================================================================== - function string_to_real( string ) - ! STRING_TO_REAL converts an arbitrary string to a real(8) - ! - ! Arguments: - ! string = character(*) containing number to convert + function str_to_real( string ) character(*), intent(in) :: string - real(8) :: string_to_real + real(8) :: str_to_real integer :: index_decimal ! index of decimal point integer :: index_exponent ! index of exponent character @@ -190,26 +204,12 @@ contains write( fmt, '("(E",I2,".",I2,")")') w, d ! Read string - read( string, fmt, iostat=readError ) string_to_real + read( string, fmt, iostat=readError ) str_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 function str_to_real end module string diff --git a/src/types.f90 b/src/types.f90 index b760cec1f..f553ad7d7 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -2,46 +2,8 @@ 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 :: uid integer :: type real(8), allocatable :: coeffs(:) integer, allocatable :: neighbor_pos(:) @@ -49,30 +11,13 @@ module types 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 :: uid 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 @@ -80,5 +25,28 @@ module types real(8) :: wgt ! particle weight end type Neutron + type :: Bank + integer :: uid ! Unique ID + real(8) :: xyz(3) ! Location of bank particle + end type Bank + + type :: Isotope + integer :: uid ! unique identifier + integer :: zaid ! ZAID, e.g. 92235 + character(3) :: xs ! cross section identifier, e.g. 70c + real(8) :: density ! density in atom/b-cm + end type Isotope + + type :: Material + integer :: uid + integer :: n_isotopes + type(Isotope), allocatable :: isotopes(:) + 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 end module types