From 1e5fb2e428865311b9b14586000fc9e9a0d3da2a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 29 Sep 2011 13:56:19 -0400 Subject: [PATCH] Added reading of meshes for tallies. Currently no logic is implemented to determine mesh bins when scoring tallies. --- src/DEPENDENCIES | 3 + src/OBJECTS | 1 + src/global.f90 | 28 ++++-- src/initialize.f90 | 17 ++++ src/input_xml.f90 | 112 +++++++++++++++++++++-- src/mesh_header.f90 | 19 ++++ src/output.f90 | 23 +++-- src/xml_tallies.f90 | 216 ++++++++++++++++++++++++++++++++++++++++++++ 8 files changed, 391 insertions(+), 28 deletions(-) create mode 100644 src/mesh_header.f90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 14709b6eab..3d11f65673 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -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 diff --git a/src/OBJECTS b/src/OBJECTS index 38d3cfdab8..1dc1cd67a3 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -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 \ diff --git a/src/global.f90 b/src/global.f90 index 96f253a30e..f3c74cdf13 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -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 diff --git a/src/initialize.f90 b/src/initialize.f90 index 1815ca2c3b..c9bc748fdb 100644 --- a/src/initialize.f90 +++ b/src/initialize.f90 @@ -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 diff --git a/src/input_xml.f90 b/src/input_xml.f90 index 4b130d1ee4..01ec7a9d29 100644 --- a/src/input_xml.f90 +++ b/src/input_xml.f90 @@ -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 must be the same as the " // & + "number of entries on ." + 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 must be the same as the " // & + "number of entries on ." + 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 diff --git a/src/mesh_header.f90 b/src/mesh_header.f90 new file mode 100644 index 0000000000..f49a5c50ab --- /dev/null +++ b/src/mesh_header.f90 @@ -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 diff --git a/src/output.f90 b/src/output.f90 index 14d040ecdc..b0a1f0843b 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -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 diff --git a/src/xml_tallies.f90 b/src/xml_tallies.f90 index b4e67a7323..3b3e9f9ff4 100644 --- a/src/xml_tallies.f90 +++ b/src/xml_tallies.f90 @@ -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)), & + '' +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)') & '' + 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)') '' call xml_close(info)