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

Added reading of meshes for tallies. Currently no logic is implemented to determine mesh bins when scoring tallies.

Browse source
This commit is contained in:
Paul Romano 2011-09-29 13:56:19 -04:00
parent e43b26b7b0
commit 1e5fb2e428
8 changed files with 391 additions and 28 deletions
Download patch file Download diff file Expand all files Collapse all files

3
src/DEPENDENCIES
View file

@ -50,6 +50,7 @@ global.o: cross_section_header.o
global.o: datatypes_header.o
global.o: geometry_header.o
global.o: material_header.o
global.o: mesh_header.o
global.o: particle_header.o
global.o: source_header.o
global.o: tally_header.o
@ -79,6 +80,7 @@ input_xml.o: datatypes.o
input_xml.o: error.o
input_xml.o: geometry_header.o
input_xml.o: global.o
input_xml.o: mesh_header.o
input_xml.o: output.o
input_xml.o: string.o
input_xml.o: tally_header.o
@ -119,6 +121,7 @@ output.o: datatypes.o
output.o: endf.o
output.o: geometry_header.o
output.o: global.o
output.o: mesh_header.o
output.o: string.o
output.o: tally_header.o

1
src/OBJECTS
View file

@ -20,6 +20,7 @@ logging.o \
main.o \
material_header.o \
mcnp_random.o \
mesh_header.o \
mpi_routines.o \
output.o \
particle_header.o \

28
src/global.f90
View file

@ -7,6 +7,7 @@ module global
use datatypes_header, only: DictionaryII, DictionaryCI
use geometry_header, only: Cell, Universe, Lattice, Surface
use material_header, only: Material
use mesh_header, only: StructuredMesh
use particle_header, only: Particle
use source_header, only: ExtSource
use tally_header, only: TallyObject, TallyMap
@ -19,29 +20,36 @@ module global
implicit none
save
! Main arrays for cells, surfaces, materials
type(Cell), allocatable, target :: cells(:)
type(Universe), allocatable, target :: universes(:)
type(Lattice), allocatable, target :: lattices(:)
type(Surface), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
type(xsData), allocatable, target :: xsdatas(:)
type(TallyObject), allocatable, target :: tallies(:)
type(TallyObject), allocatable, target :: tallies_global(:)
! Main arrays
type(Cell), allocatable, target :: cells(:)
type(Universe), allocatable, target :: universes(:)
type(Lattice), allocatable, target :: lattices(:)
type(Surface), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
type(xsData), allocatable, target :: xsdatas(:)
type(StructuredMesh), allocatable, target :: meshes(:)
type(TallyObject), allocatable, target :: tallies(:)
type(TallyObject), allocatable, target :: tallies_global(:)
! Size of main arrays
integer :: n_cells ! # of cells
integer :: n_universes ! # of universes
integer :: n_lattices ! # of lattices
integer :: n_surfaces ! # of surfaces
integer :: n_materials ! # of materials
integer :: n_meshes ! # of structured meshes
integer :: n_tallies ! # of tallies
integer :: n_tallies_global ! # of global tallies
! These dictionaries provide a fast lookup mechanism
! These dictionaries provide a fast lookup mechanism -- the key is the
! user-specified identifier and the value is the index in the corresponding
! array
type(DictionaryII), pointer :: cell_dict
type(DictionaryII), pointer :: universe_dict
type(DictionaryII), pointer :: lattice_dict
type(DictionaryII), pointer :: surface_dict
type(DictionaryII), pointer :: material_dict
type(DictionaryII), pointer :: mesh_dict
type(DictionaryII), pointer :: tally_dict
type(DictionaryCI), pointer :: xsdata_dict
type(DictionaryCI), pointer :: nuclide_dict

17
src/initialize.f90
View file

@ -160,6 +160,7 @@ contains
call dict_create(lattice_dict)
call dict_create(surface_dict)
call dict_create(material_dict)
call dict_create(mesh_dict)
call dict_create(tally_dict)
call dict_create(cells_in_univ_dict)
@ -411,6 +412,22 @@ contains
end if
end do
end if
! =======================================================================
! ADJUST MESH INDICES FOR EACH TALLY
if (t % n_bins(T_MESH) > 0) then
do j = 1, size(t % mesh_bins)
uid = t % mesh_bins(j) % scalar
if (dict_has_key(mesh_dict, uid)) then
t % mesh_bins(j) % scalar = dict_get_key(mesh_dict, uid)
else
msg = "Could not find mesh " // trim(int_to_str(uid)) // &
& " specified on tally " // trim(int_to_str(t % uid))
call fatal_error(msg)
end if
end do
end if
end do
end subroutine adjust_indices

112
src/input_xml.f90
View file

@ -3,9 +3,10 @@ module input_xml
use constants
use datatypes, only: dict_create, dict_add_key, dict_has_key, &
dict_get_key
use error, only: fatal_error
use error, only: fatal_error, warning
use geometry_header, only: Cell, Surface, Lattice
use global
use mesh_header, only: StructuredMesh
use output, only: message
use string, only: lower_case, int_to_str, str_to_int, str_to_real, &
split_string
@ -501,15 +502,19 @@ contains
use xml_data_tallies_t
integer :: i ! loop over user-specified tallies
integer :: j ! loop over words
integer :: n_words
integer :: i ! loop over user-specified tallies
integer :: j ! loop over words
integer :: uid ! user-specified identifier
integer :: index ! index in meshes array
integer :: n ! size of arrays in mesh specification
integer :: n_words ! number of words read
logical :: file_exists ! does tallies.xml file exist?
character(MAX_LINE_LEN) :: filename
character(MAX_LINE_LEN) :: msg
character(MAX_WORD_LEN) :: word
character(MAX_WORD_LEN) :: words(MAX_WORDS)
logical :: file_exists
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
! Check if tallies.xml exists
filename = trim(path_input) // "tallies.xml"
@ -527,9 +532,87 @@ contains
! Parse tallies.xml file
call read_xml_file_tallies_t(filename)
! Allocate cells array
n_tallies = size(tally_)
allocate(tallies(n_tallies))
! ==========================================================================
! DETERMINE SIZE OF ARRAYS AND ALLOCATE
! Check for meshes and allocate
if (.not. associated(mesh_)) then
n_meshes = 0
else
n_meshes = size(mesh_)
allocate(meshes(n_meshes))
end if
! Allocate tallies array
if (.not. associated(tally_)) then
n_tallies = 0
msg = "No tallies present in tallies.xml file!"
call warning(msg)
else
n_tallies = size(tally_)
allocate(tallies(n_tallies))
end if
! ==========================================================================
! READ MESH DATA
do i = 1, n_meshes
m => meshes(i)
! copy mesh uid
m % uid = mesh_(i) % id
! Read mesh type
word = mesh_(i) % type
call lower_case(word)
select case (trim(word))
case ('rect', 'rectangle', 'rectangular')
m % type = LATTICE_RECT
case ('hex', 'hexagon', 'hexagonal')
m % type = LATTICE_HEX
case default
msg = "Invalid mesh type: " // trim(mesh_(i) % type)
call fatal_error(msg)
end select
! Determine number of dimensions for mesh
n = size(mesh_(i) % dimension)
if (n /= 2 .and. n /= 3) then
msg = "Mesh must be two or three dimensions."
call fatal_error(msg)
end if
m % n_dimension = n
! Allocate attribute arrays
allocate(m % dimension(n))
allocate(m % origin(n))
allocate(m % width(n))
! Read dimensions in each direction
m % dimension = mesh_(i) % dimension
! Read mesh origin location
if (m % n_dimension /= size(mesh_(i) % origin)) then
msg = "Number of entries on <origin> must be the same as the " // &
"number of entries on <dimension>."
call fatal_error(msg)
end if
m % origin = mesh_(i) % origin
! Read mesh widths
if (size(mesh_(i) % width) /= size(mesh_(i) % origin)) then
msg = "Number of entries on <width> must be the same as the " // &
"number of entries on <origin>."
call fatal_error(msg)
end if
m % width = mesh_(i) % width
! Add mesh to dictionary
call dict_add_key(mesh_dict, m % uid, i)
end do
! ==========================================================================
! READ TALLY DATA
do i = 1, n_tallies
t => tallies(i)
@ -601,6 +684,19 @@ contains
allocate(t % mesh_bins(n_words))
do j = 1, n_words
t % mesh_bins(j) % scalar = str_to_int(words(j))
! Determine index in mesh array for this bin
uid = t % mesh_bins(j) % scalar
if (dict_has_key(mesh_dict, uid)) then
index = dict_get_key(mesh_dict, uid)
m => meshes(index)
else
msg = "Could not find mesh " // trim(int_to_str(uid)) // &
" specified on tally " // trim(int_to_str(t % uid))
call fatal_error(msg)
end if
t % n_bins(T_MESH) = t % n_bins(T_MESH) + product(m % dimension)
end do
end if

19
src/mesh_header.f90 Normal file
View file

@ -0,0 +1,19 @@
module mesh_header
implicit none
!===============================================================================
! STRUCTUREDMESH represents a tesslation of n-dimensional Euclidean space by
! congruent squares or cubes
!===============================================================================
type StructuredMesh
integer :: uid
integer :: type
integer :: n_dimension
integer, allocatable :: dimension(:)
integer, allocatable :: origin(:)
integer, allocatable :: width(:)
end type StructuredMesh
end module mesh_header

23
src/output.f90
View file

@ -8,6 +8,7 @@ module output
use endf, only: reaction_name
use geometry_header, only: Cell, Universe, Surface
use global
use mesh_header, only: StructuredMesh
use string, only: upper_case, int_to_str, real_to_str
use tally_header, only: TallyObject
@ -506,13 +507,14 @@ contains
type(TallyObject), pointer :: t
integer :: i
integer :: uid
character(MAX_LINE_LEN) :: string
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
type(Material), pointer :: m => null()
integer :: i
integer :: uid
character(MAX_LINE_LEN) :: string
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
type(Material), pointer :: m => null()
type(StructuredMesh), pointer :: mesh => null()
write(ou,*) 'Tally ' // int_to_str(t % uid)
@ -556,11 +558,12 @@ contains
write(ou, *) ' Material Bins:' // trim(string)
end if
if (associated(t % mesh_bins)) then
if (t % n_bins(T_MESH) > 0) then
string = ""
do i = 1, size(t % mesh_bins)
string = trim(string) // ' ' // trim(int_to_str(&
t % mesh_bins(i) % scalar))
uid = t % mesh_bins(i) % scalar
mesh => meshes(uid)
string = trim(string) // ' ' // trim(int_to_str(mesh % uid))
end do
write(ou, *) ' Mesh Bins:' // trim(string)
end if

216
src/xml_tallies.f90
View file

@ -6,6 +6,14 @@ module xml_data_tallies_t
integer, private :: lurep_
logical, private :: strict_
type mesh_xml
integer :: id
character(len=12) :: type
integer, dimension(:), pointer :: dimension => null()
real(kind=kind(1.0d0)), dimension(:), pointer :: origin => null()
real(kind=kind(1.0d0)), dimension(:), pointer :: width => null()
end type mesh_xml
type filter_xml
character(len=250) :: cell
character(len=250) :: surface
@ -24,8 +32,204 @@ type tally_xml
character(len=250) :: reactions
character(len=250) :: nuclides
end type tally_xml
type(mesh_xml), dimension(:), pointer :: mesh_ => null()
type(tally_xml), dimension(:), pointer :: tally_ => null()
contains
subroutine read_xml_type_mesh_xml_array( &
info, tag, endtag, attribs, noattribs, data, nodata, &
dvar, has_dvar )
type(XML_PARSE) :: info
character(len=*), intent(inout) :: tag
logical, intent(inout) :: endtag
character(len=*), dimension(:,:), intent(inout) :: attribs
integer, intent(inout) :: noattribs
character(len=*), dimension(:), intent(inout) :: data
integer, intent(inout) :: nodata
type(mesh_xml), dimension(:), pointer :: dvar
logical, intent(inout) :: has_dvar
integer :: newsize
type(mesh_xml), dimension(:), pointer :: newvar
newsize = size(dvar) + 1
allocate( newvar(1:newsize) )
newvar(1:newsize-1) = dvar
deallocate( dvar )
dvar => newvar
call read_xml_type_mesh_xml( info, tag, endtag, attribs, noattribs, data, nodata, &
dvar(newsize), has_dvar )
end subroutine read_xml_type_mesh_xml_array
subroutine read_xml_type_mesh_xml( info, starttag, endtag, attribs, noattribs, data, nodata, &
dvar, has_dvar )
type(XML_PARSE) :: info
character(len=*), intent(in) :: starttag
logical, intent(inout) :: endtag
character(len=*), dimension(:,:), intent(inout) :: attribs
integer, intent(inout) :: noattribs
character(len=*), dimension(:), intent(inout) :: data
integer, intent(inout) :: nodata
type(mesh_xml), intent(inout) :: dvar
logical, intent(inout) :: has_dvar
integer :: att_
integer :: noatt_
logical :: error
logical :: endtag_org
character(len=len(starttag)) :: tag
logical :: has_id
logical :: has_type
logical :: has_dimension
logical :: has_origin
logical :: has_width
has_id = .false.
has_type = .false.
has_dimension = .false.
has_origin = .false.
has_width = .false.
call init_xml_type_mesh_xml(dvar)
has_dvar = .true.
error = .false.
att_ = 0
noatt_ = noattribs+1
endtag_org = endtag
do
if ( nodata /= 0 ) then
noattribs = 0
tag = starttag
elseif ( att_ < noatt_ .and. noatt_ > 1 ) then
att_ = att_ + 1
if ( att_ <= noatt_-1 ) then
tag = attribs(1,att_)
data(1) = attribs(2,att_)
noattribs = 0
nodata = 1
endtag = .false.
else
tag = starttag
noattribs = 0
nodata = 0
endtag = .true.
cycle
endif
else
if ( endtag_org ) then
return
else
call xml_get( info, tag, endtag, attribs, noattribs, data, nodata )
if ( xml_error(info) ) then
write(lurep_,*) 'Error reading input file!'
error = .true.
return
endif
endif
endif
if ( endtag .and. tag == starttag ) then
exit
endif
if ( endtag .and. noattribs == 0 ) then
if ( xml_ok(info) ) then
cycle
else
exit
endif
endif
select case( tag )
case('id')
call read_xml_integer( &
info, tag, endtag, attribs, noattribs, data, nodata, &
dvar%id, has_id )
case('type')
call read_xml_word( &
info, tag, endtag, attribs, noattribs, data, nodata, &
dvar%type, has_type )
case('dimension')
call read_xml_integer_array( &
info, tag, endtag, attribs, noattribs, data, nodata, &
dvar%dimension, has_dimension )
case('origin')
call read_xml_double_array( &
info, tag, endtag, attribs, noattribs, data, nodata, &
dvar%origin, has_origin )
case('width')
call read_xml_double_array( &
info, tag, endtag, attribs, noattribs, data, nodata, &
dvar%width, has_width )
case ('comment', '!--')
! Simply ignore
case default
if ( strict_ ) then
error = .true.
call xml_report_errors( info, &
'Unknown or wrongly placed tag: ' // trim(tag))
endif
end select
nodata = 0
if ( .not. xml_ok(info) ) exit
end do
if ( .not. has_id ) then
has_dvar = .false.
call xml_report_errors(info, 'Missing data on id')
endif
if ( .not. has_type ) then
has_dvar = .false.
call xml_report_errors(info, 'Missing data on type')
endif
if ( .not. has_dimension ) then
has_dvar = .false.
call xml_report_errors(info, 'Missing data on dimension')
endif
if ( .not. has_origin ) then
has_dvar = .false.
call xml_report_errors(info, 'Missing data on origin')
endif
if ( .not. has_width ) then
has_dvar = .false.
call xml_report_errors(info, 'Missing data on width')
endif
end subroutine read_xml_type_mesh_xml
subroutine init_xml_type_mesh_xml_array( dvar )
type(mesh_xml), dimension(:), pointer :: dvar
if ( associated( dvar ) ) then
deallocate( dvar )
endif
allocate( dvar(0) )
end subroutine init_xml_type_mesh_xml_array
subroutine init_xml_type_mesh_xml(dvar)
type(mesh_xml) :: dvar
end subroutine init_xml_type_mesh_xml
subroutine write_xml_type_mesh_xml_array( &
info, tag, indent, dvar )
type(XML_PARSE) :: info
character(len=*), intent(in) :: tag
integer :: indent
type(mesh_xml), dimension(:) :: dvar
integer :: i
do i = 1,size(dvar)
call write_xml_type_mesh_xml( info, tag, indent, dvar(i) )
enddo
end subroutine write_xml_type_mesh_xml_array
subroutine write_xml_type_mesh_xml( &
info, tag, indent, dvar )
type(XML_PARSE) :: info
character(len=*), intent(in) :: tag
integer :: indent
type(mesh_xml) :: dvar
character(len=100) :: indentation
indentation = ' '
write(info%lun, '(4a)' ) indentation(1:min(indent,100)),&
'<',trim(tag), '>'
call write_to_xml_integer( info, 'id', indent+3, dvar%id)
call write_to_xml_word( info, 'type', indent+3, dvar%type)
call write_to_xml_integer_array( info, 'dimension', indent+3, dvar%dimension)
call write_to_xml_double_array( info, 'origin', indent+3, dvar%origin)
call write_to_xml_double_array( info, 'width', indent+3, dvar%width)
write(info%lun,'(4a)') indentation(1:min(indent,100)), &
'</' //trim(tag) // '>'
end subroutine write_xml_type_mesh_xml
subroutine read_xml_type_filter_xml_array( &
info, tag, endtag, attribs, noattribs, data, nodata, &
dvar, has_dvar )
@ -430,7 +634,10 @@ subroutine read_xml_file_tallies_t(fname, lurep, errout)
integer :: noattribs
character(len=200), dimension(1:100) :: data
integer :: nodata
logical :: has_mesh_
logical :: has_tally_
has_mesh_ = .false.
allocate(mesh_(0))
has_tally_ = .false.
allocate(tally_(0))
@ -474,6 +681,10 @@ subroutine read_xml_file_tallies_t(fname, lurep, errout)
endif
endif
select case( tag )
case('mesh')
call read_xml_type_mesh_xml_array( &
info, tag, endtag, attribs, noattribs, data, nodata, &
mesh_, has_mesh_ )
case('tally')
call read_xml_type_tally_xml_array( &
info, tag, endtag, attribs, noattribs, data, nodata, &
@ -490,6 +701,10 @@ subroutine read_xml_file_tallies_t(fname, lurep, errout)
nodata = 0
if ( .not. xml_ok(info) ) exit
end do
if ( .not. has_mesh_ ) then
error = .true.
call xml_report_errors(info, 'Missing data on mesh_')
endif
if ( .not. has_tally_ ) then
error = .true.
call xml_report_errors(info, 'Missing data on tally_')
@ -511,6 +726,7 @@ subroutine write_xml_file_tallies_t(fname, lurep)
endif
write(info%lun,'(a)') &
'<tallies>'
call write_xml_type_mesh_xml_array( info, 'mesh', indent+3, mesh_)
call write_xml_type_tally_xml_array( info, 'tally', indent+3, tally_)
write(info%lun,'(a)') '</tallies>'
call xml_close(info)