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

Changed allocation of main arrays, many formatting and structural changes.

Browse source
This commit is contained in:
Paul Romano 2011-01-23 05:05:58 +00:00
parent 4d044c6a61
commit 6d5cafbcb7
12 changed files with 1272 additions and 316 deletions
Download patch file Download diff file Expand all files Collapse all files

19
ChangeLog Normal file
View file

@ -0,0 +1,19 @@
2011-01-22 Paul Romano <romano7@mit.edu>
* 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

8
src/Makefile
View file

@ -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

377
src/fileio.f90
View file

@ -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

33
src/geometry.f90
View file

@ -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 )

81
src/global.f90
View file

@ -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

14
src/input_sample Normal file
View file

@ -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

77
src/main.f90
View file

@ -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

614
src/mcnp_random.f90 Normal file
View file

@ -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

146
src/output.f90
View file

@ -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

33
src/source.f90 Normal file
View file

@ -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

104
src/string.f90
View file

@ -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

82
src/types.f90
View file

@ -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