diff --git a/docs/source/developers.rst b/docs/source/developers.rst index d4e9570a8d..56ce5e14ad 100644 --- a/docs/source/developers.rst +++ b/docs/source/developers.rst @@ -12,6 +12,7 @@ Active development of the OpenMC Monte Carlo code is currently led by: * `Adam Nelson `_ * `Jon Walsh `_ * `Sterling Harper `_ +* `Will Boyd `_ * `Benoit Forget `_ * `Kord Smith `_ * `Andrew Siegel `_ diff --git a/docs/source/quickinstall.rst b/docs/source/quickinstall.rst index 3fb9036d39..38beb6cca5 100644 --- a/docs/source/quickinstall.rst +++ b/docs/source/quickinstall.rst @@ -42,18 +42,21 @@ the following commands in a terminal: .. code-block:: sh git clone git://github.com/mit-crpg/openmc.git - cd openmc/src + cd openmc git checkout -b master origin/master + mkdir src/build + cd src/build + cmake .. make sudo make install This will build an executable named ``openmc`` and install it (by default in -/usr/local/bin). If you do not have administrator privileges, the last command -can be replaced with a local install, e.g. +/usr/local/bin). If you do not have administrator privileges, the cmake command +should specify an installation directory where you have write access, e.g. .. code-block:: sh - make install -e prefix=$HOME/.local + cmake -DCMAKE_INSTALL_PREFIX=$HOME/.local .. .. _GitHub: https://github.com/mit-crpg/openmc .. _git: http://git-scm.com diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index c1231e9a8e..910cbb0503 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -708,6 +708,12 @@ Each ```` element can have the following attributes or sub-elements: *Default*: None + :name: + An optional string name to identify the surface in summary output + files. This string is limited to 52 characters for formatting purposes. + + *Default*: "" + :type: The type of the surfaces. This can be "x-plane", "y-plane", "z-plane", "plane", "x-cylinder", "y-cylinder", "z-cylinder", or "sphere". @@ -788,6 +794,12 @@ Each ```` element can have the following attributes or sub-elements: *Default*: None + :name: + An optional string name to identify the cell in summary output files. + This string is limmited to 52 characters for formatting purposes. + + *Default*: "" + :universe: The ``id`` of the universe that this cell is contained in. @@ -853,6 +865,12 @@ the following attributes or sub-elements: :id: A unique integer that can be used to identify the lattice. + :name: + An optional string name to identify the lattice in summary output + files. This string is limited to 52 characters for formatting purposes. + + *Default*: "" + :dimension: Two or three integers representing the number of lattice cells in the x- and y- (and z-) directions, respectively. @@ -912,6 +930,12 @@ the following attributes or sub-elements: :id: A unique integer that can be used to identify the lattice. + :name: + An optional string name to identify the hex_lattice in summary output + files. This string is limited to 52 characters for formatting purposes. + + *Default*: "" + :n_rings: An integer representing the number of radial ring positions in the xy-plane. Note that this number includes the degenerate center ring which only has one @@ -989,6 +1013,12 @@ Each ``material`` element can have the following attributes or sub-elements: :id: A unique integer that can be used to identify the material. + :name: + An optional string name to identify the material in summary output + files. This string is limited to 52 characters for formatting purposes. + + *Default*: "" + :density: An element with attributes/sub-elements called ``value`` and ``units``. The ``value`` attribute is the numeric value of the density while the ``units`` @@ -1100,10 +1130,11 @@ and ````. The ```` element accepts the following sub-elements: - :label: - This is an optional sub-element specifying the name of this tally to be used - for output purposes. This string is limited to 52 characters for formatting - purposes. + :name: + An optional string name to identify the tally in summary output + files. This string is limited to 52 characters for formatting purposes. + + *Default*: "" :filter: Specify a filter that restricts contributions to the tally to particles diff --git a/examples/python/pincell/build-xml.py b/examples/python/pincell/build-xml.py index c79fec7fe9..9338aff0e9 100644 --- a/examples/python/pincell/build-xml.py +++ b/examples/python/pincell/build-xml.py @@ -200,7 +200,7 @@ mesh_filter = openmc.Filter() mesh_filter.mesh = mesh # Instantiate the Tally -tally = openmc.Tally(tally_id=1) +tally = openmc.Tally(tally_id=1, name='tally 1') tally.add_filter(energy_filter) tally.add_filter(mesh_filter) tally.add_score('flux') diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d3785aa6a3..5fb11848db 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -261,15 +261,12 @@ set_target_properties(${program} PROPERTIES #=============================================================================== install(TARGETS ${program} RUNTIME DESTINATION bin) -install(PROGRAMS utils/statepoint_cmp.py +install(PROGRAMS utils/plot_mesh_tally.py DESTINATION bin - RENAME statepoint_cmp) + RENAME plot_mesh_tally) install(PROGRAMS utils/statepoint_histogram.py DESTINATION bin RENAME statepoint_histogram) -install(PROGRAMS utils/statepoint_meshplot.py - DESTINATION bin - RENAME statepoint_meshplot) install(PROGRAMS utils/update_inputs.py DESTINATION bin RENAME update_inputs) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 578ecb9678..11059c48f2 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -460,8 +460,8 @@ contains if (i == 1) then - ! Set label - t % label = "CMFD flux, total, scatter-1" + ! Set name + t % name = "CMFD flux, total, scatter-1" ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG @@ -491,8 +491,8 @@ contains else if (i == 2) then - ! Set label - t % label = "CMFD neutron production" + ! Set name + t % name = "CMFD neutron production" ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG @@ -536,8 +536,8 @@ contains else if (i == 3) then - ! Set label - t % label = "CMFD surface currents" + ! Set name + t % name = "CMFD surface currents" ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 6f6cfaecec..b6cf43119e 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -22,6 +22,7 @@ module geometry_header type, abstract :: Lattice integer :: id ! Universe number for lattice + character(len=52) :: name = "" ! User-defined name real(8), allocatable :: pitch(:) ! Pitch along each axis integer, allocatable :: universes(:,:,:) ! Specified universes integer :: outside ! Material to fill area outside @@ -116,6 +117,7 @@ module geometry_header type Surface integer :: id ! Unique ID + character(len=52) :: name = "" ! User-defined name integer :: type ! Type of surface real(8), allocatable :: coeffs(:) ! Definition of surface integer, allocatable :: & @@ -130,6 +132,7 @@ module geometry_header type Cell integer :: id ! Unique ID + character(len=52) :: name = "" ! User-defined name integer :: type ! Type of cell (normal, universe, lattice) integer :: universe ! universe # this cell is in integer :: fill ! universe # filling this cell diff --git a/src/hdf5_summary.F90 b/src/hdf5_summary.F90 index c4ad25eb5b..4b203501be 100644 --- a/src/hdf5_summary.F90 +++ b/src/hdf5_summary.F90 @@ -127,6 +127,10 @@ contains call su % write_data(i, "index", & group="geometry/cells/cell " // trim(to_str(c % id))) + ! Write name for this cell + call su % write_data(c % name, "name", & + group="geometry/cells/cell " // trim(to_str(c % id))) + ! Write universe for this cell call su % write_data(universes(c % universe) % id, "universe", & group="geometry/cells/cell " // trim(to_str(c % id))) @@ -200,6 +204,10 @@ contains call su % write_data(i, "index", & group="geometry/surfaces/surface " // trim(to_str(s % id))) + ! Write name for this surface + call su % write_data(s % name, "name", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + ! Write surface type select case (s % type) case (SURF_PX) @@ -311,6 +319,10 @@ contains call su % write_data(i, "index", & group="geometry/lattices/lattice " // trim(to_str(lat % id))) + ! Write name for this lattice + call su % write_data(lat % name, "name", & + group="geometry/lattices/lattice " // trim(to_str(lat % id))) + ! Write lattice type select type (lat) type is (RectLattice) @@ -437,6 +449,10 @@ contains call su % write_data(i, "index", & group="materials/material " // trim(to_str(m % id))) + ! Write name for this material + call su % write_data(m % name, "name", & + group="materials/material " // trim(to_str(m % id))) + ! Write atom density with units call su % write_data(m % density, "atom_density", & group="materials/material " // trim(to_str(m % id))) @@ -531,6 +547,14 @@ contains ! Get pointer to tally t => tallies(i) + ! Write the name for this tally + call su % write_data(len(t % name), "name_size", & + group="tallies/tally " // trim(to_str(t % id))) + if (len(t % name) > 0) then + call su % write_data(t % name, "name", & + group="tallies/tally " // trim(to_str(t % id))) + endif + ! Write size of each tally call su % write_data(t % total_score_bins, "total_score_bins", & group="tallies/tally " // trim(to_str(t % id))) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 01379a231c..93129d0544 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -975,6 +975,12 @@ contains else call fatal_error("Must specify id of cell in geometry XML file.") end if + + ! Copy cell name + if (check_for_node(node_cell, "name")) then + call get_node_value(node_cell, "name", c % name) + end if + if (check_for_node(node_cell, "universe")) then call get_node_value(node_cell, "universe", c % universe) else @@ -1154,6 +1160,11 @@ contains &// to_str(s % id)) end if + ! Copy surface name + if (check_for_node(node_surf, "name")) then + call get_node_value(node_surf, "name", s % name) + end if + ! Copy and interpret surface type word = '' if (check_for_node(node_surf, "type")) & @@ -1276,6 +1287,11 @@ contains &// to_str(lat % id)) end if + ! Copy lattice name + if (check_for_node(node_lat, "name")) then + call get_node_value(node_lat, "name", lat % name) + end if + ! Read number of lattice cells in each dimension n = get_arraysize_integer(node_lat, "dimension") if (n == 2) then @@ -1399,6 +1415,11 @@ contains &// to_str(lat % id)) end if + ! Copy lattice name + if (check_for_node(node_lat, "name")) then + call get_node_value(node_lat, "name", lat % name) + end if + ! Read number of lattice cells in each dimension call get_node_value(node_lat, "n_rings", lat % n_rings) if (check_for_node(node_lat, "n_axial")) then @@ -1645,6 +1666,11 @@ contains &// to_str(mat % id)) end if + ! Copy material name + if (check_for_node(node_mat, "name")) then + call get_node_value(node_mat, "name", mat % name) + end if + if (run_mode == MODE_PLOTTING) then ! add to the dictionary and skip xs processing call material_dict % add_key(mat % id, i) @@ -2232,10 +2258,9 @@ contains &// to_str(t % id)) end if - ! Copy tally label - t % label = '' - if (check_for_node(node_tal, "label")) & - call get_node_value(node_tal, "label", t % label) + ! Copy tally name + if (check_for_node(node_tal, "name")) & + call get_node_value(node_tal, "name", t % name) ! ======================================================================= ! READ DATA FOR FILTERS diff --git a/src/material_header.F90 b/src/material_header.F90 index b992ed16ff..e36250176c 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -8,6 +8,7 @@ module material_header type Material integer :: id ! unique identifier + character(len=52) :: name = "" ! User-defined name integer :: n_nuclides ! number of nuclides integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm diff --git a/src/mpiio_interface.F90 b/src/mpiio_interface.F90 index 97ca73fb3c..30093adec9 100644 --- a/src/mpiio_interface.F90 +++ b/src/mpiio_interface.F90 @@ -21,7 +21,6 @@ module mpiio_interface module procedure mpi_write_integer_4Darray module procedure mpi_write_long module procedure mpi_write_string - !module procedure mpi_write_string_1Darray end interface mpi_write_data ! Generic HDF5 read procedure interface @@ -38,7 +37,6 @@ module mpiio_interface module procedure mpi_read_integer_4Darray module procedure mpi_read_long module procedure mpi_read_string - !module procedure mpi_read_string_1Darray end interface mpi_read_data contains diff --git a/src/output.F90 b/src/output.F90 index bd86835d10..3cea71de6a 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -369,6 +369,9 @@ contains ! Write user-specified id for cell write(unit_,*) 'Cell ' // to_str(c % id) + ! Write user-specified name for cell + write(unit_,*) ' Name = ' // c % name + ! Find index in cells array and write index_cell = cell_dict % get_key(c % id) write(unit_,*) ' Array Index = ' // to_str(index_cell) @@ -487,6 +490,9 @@ contains ! Write information about lattice write(unit_,*) 'Lattice ' // to_str(lat % id) + ! Write user-specified name for lattice + write(unit_,*) ' Name = ' // lat % name + select type(lat) type is (RectLattice) ! Write dimension of lattice. @@ -572,6 +578,9 @@ contains ! Write user-specified id of surface write(unit_,*) 'Surface ' // to_str(surf % id) + ! Write user-specified name for surface + write(unit_,*) ' Name = ' // surf % name + ! Write type of surface select case (surf % type) case (SURF_PX) @@ -668,6 +677,9 @@ contains ! Write identifier for material write(unit_,*) 'Material ' // to_str(mat % id) + ! Write user-specified name for material + write(unit_,*) ' Name = ' // mat % name + ! Write total atom density in atom/b-cm write(unit_,*) ' Atom Density = ' // trim(to_str(mat % density)) & // ' atom/b-cm' @@ -1767,12 +1779,12 @@ contains end if ! Write header block - if (t % label == "") then + if (t % name == "") then call header("TALLY " // trim(to_str(t % id)), unit=UNIT_TALLY, & level=3) else call header("TALLY " // trim(to_str(t % id)) // ": " & - // trim(t % label), unit=UNIT_TALLY, level=3) + // trim(t % name), unit=UNIT_TALLY, level=3) endif ! Handle surface current tallies separately diff --git a/src/output_interface.F90 b/src/output_interface.F90 index 6d43de21f7..d8134aa583 100644 --- a/src/output_interface.F90 +++ b/src/output_interface.F90 @@ -1761,91 +1761,6 @@ contains #endif end subroutine read_string -!=============================================================================== -! WRITE_STRING_1DARRAY writes 1-D string data -!=============================================================================== - - subroutine write_string_1Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length ! length of array to write - character(*), intent(in) :: buffer(:) ! data to write - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - integer :: n - logical :: collect_ - - ! Get string length - n = len_trim(buffer(1)) - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - -#ifdef HDF5 - ! Check if HDF5 group should be created/opened - if (present(group)) then - call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) - else - self % hdf5_grp = self % hdf5_fh - endif - if (present(group)) call hdf5_close_group(self % hdf5_grp) -#endif - - end subroutine write_string_1Darray - -!=============================================================================== -! READ_STRING_1DARRAY reads 1-D string data -!=============================================================================== - - subroutine read_string_1Darray(self, buffer, name, group, length, collect) - - integer, intent(in) :: length ! length of array to write - character(*), intent(inout) :: buffer(:) ! data to write - character(*), intent(in) :: name ! name of data - character(*), intent(in), optional :: group ! HDF5 group name - logical, intent(in), optional :: collect ! collective I/O - class(BinaryOutput) :: self - - character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name - character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name - integer :: n - logical :: collect_ - - ! Get string length - n = len(buffer(1)) - - ! Set name - name_ = trim(name) - - ! Set group - if (present(group)) then - group_ = trim(group) - end if - - ! Set up collective vs. independent I/O - if (present(collect)) then - collect_ = collect - else - collect_ = .true. - end if - - end subroutine read_string_1Darray !=============================================================================== ! WRITE_ATTRIBUTE_STRING diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc index 801c902c9e..cebbc5b6a3 100644 --- a/src/relaxng/geometry.rnc +++ b/src/relaxng/geometry.rnc @@ -1,6 +1,8 @@ element geometry { element cell { (element id { xsd:int } | attribute id { xsd:int }) & + (element name { xsd:string { maxLength="52" } } | + attribute name { xsd:string { maxLength="52" } })? & (element universe { xsd:int } | attribute universe { xsd:int })? & ( (element fill { xsd:int } | attribute fill { xsd:int }) | @@ -14,6 +16,8 @@ element geometry { & element surface { (element id { xsd:int } | attribute id { xsd:int }) & + (element name { xsd:string { maxLength="52" } } | + attribute name { xsd:string { maxLength="52" } })? & (element type { xsd:string { maxLength = "15" } } | attribute type { xsd:string { maxLength = "15" } }) & (element coeffs { list { xsd:double+ } } | attribute coeffs { list { xsd:double+ } }) & @@ -23,6 +27,8 @@ element geometry { & element lattice { (element id { xsd:int } | attribute id { xsd:int }) & + (element name { xsd:string { maxLength="52" } } | + attribute name { xsd:string { maxLength="52" } })? & (element dimension { list { xsd:positiveInteger+ } } | attribute dimension { list { xsd:positiveInteger+ } }) & (element lower_left { list { xsd:double+ } } | attribute lower_left { list { xsd:double+ } }) & @@ -33,6 +39,8 @@ element geometry { & element hex_lattice { (element id { xsd:int } | attribute id { xsd:int }) & + (element name { xsd:string { maxLength="52" } } | + attribute name { xsd:string { maxLength="52" } })? & (element n_rings { xsd:int } | attribute n_rings { xsd:int }) & (element n_axial { xsd:int } | attribute n_axial { xsd:int })? & (element center { list { xsd:double+ } } | attribute center { list { xsd:double+ } }) & diff --git a/src/relaxng/geometry.rng b/src/relaxng/geometry.rng index 9ac3c9c598..fdbf74cbd5 100644 --- a/src/relaxng/geometry.rng +++ b/src/relaxng/geometry.rng @@ -12,6 +12,20 @@ + + + + + 52 + + + + + 52 + + + + @@ -114,6 +128,20 @@ + + + + + 52 + + + + + 52 + + + + @@ -174,6 +202,20 @@ + + + + + 52 + + + + + 52 + + + + @@ -262,6 +304,20 @@ + + + + + 52 + + + + + 52 + + + + diff --git a/src/relaxng/materials.rnc b/src/relaxng/materials.rnc index d497fa745f..ae85654972 100644 --- a/src/relaxng/materials.rnc +++ b/src/relaxng/materials.rnc @@ -1,6 +1,8 @@ element materials { element material { (element id { xsd:int } | attribute id { xsd:int }) & + (element name { xsd:string { maxLength="52" } } | + attribute name { xsd:string { maxLength="52" } })? & element density { (element value { xsd:double } | attribute value { xsd:double })? & diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index 09a5a6eb85..6911dab2f8 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -17,8 +17,8 @@ element tallies { element tally { (element id { xsd:int } | attribute id { xsd:int }) & - (element label { xsd:string { maxLength="52" } } | - attribute label { xsd:string { maxLength="52" } })? & + (element name { xsd:string { maxLength="52" } } | + attribute name { xsd:string { maxLength="52" } })? & (element estimator { ( "analog" | "tracklength" ) } | attribute estimator { ( "analog" | "tracklength" ) })? & element filter { diff --git a/src/relaxng/tallies.rng b/src/relaxng/tallies.rng index 60263676be..a3668925ad 100644 --- a/src/relaxng/tallies.rng +++ b/src/relaxng/tallies.rng @@ -108,12 +108,12 @@ - + 52 - + 52 diff --git a/src/state_point.F90 b/src/state_point.F90 index e30be0ea31..ea5d1ed46b 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -236,13 +236,6 @@ contains ! Get pointer to tally tally => tallies(i) - call sp % write_data(len(tally % label), "label_size", & - group="tallies/tally " // trim(to_str(tally % id))) - if (len(tally % label) > 0) then - call sp % write_data(tally % label, "label", & - group="tallies/tally " // trim(to_str(tally % id))) - endif - call sp % write_data(tally % estimator, "estimator", & group="tallies/tally " // trim(to_str(tally % id))) call sp % write_data(tally % n_realizations, "n_realizations", & @@ -300,7 +293,7 @@ contains MOMENT_LOOP: do j = 1, tally % n_user_score_bins select case(tally % score_bins(k)) case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) - moment_name = 'P' // to_str(tally % moment_order(k)) + moment_name = 'P' // trim(to_str(tally % moment_order(k))) call sp % write_data(moment_name, "order" // trim(to_str(k)), & group="tallies/tally " // trim(to_str(tally % id)) // & "/moments") @@ -649,7 +642,7 @@ contains character(MAX_FILE_LEN) :: path_temp character(19) :: current_time - character(52) :: label + character(52) :: name integer :: i, j, k integer :: length(4) integer :: int_array(3) @@ -826,13 +819,6 @@ contains tally => tallies(i) curr_key = key_array(id_array(i)) - call sp % read_data(j, "label_size", group="tallies/tally " // & - trim(to_str(curr_key))) - if (j > 0) then - call sp % read_data(label, "label", group="tallies/tally " // & - trim(to_str(curr_key))) - end if - call sp % read_data(tally % estimator, "estimator", & group="tallies/tally " // trim(to_str(curr_key))) call sp % read_data(tally % n_realizations, "n_realizations", & diff --git a/src/tally_header.F90 b/src/tally_header.F90 index e440f472e8..1487f1ddd7 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -71,7 +71,7 @@ module tally_header ! Basic data integer :: id ! user-defined identifier - character(len=52) :: label = "" ! user-defined label + character(len=52) :: name = "" ! user-defined name integer :: type ! volume, surface current integer :: estimator ! collision, track-length real(8) :: volume ! volume of region @@ -159,7 +159,7 @@ module tally_header ! This routine will go through each item in TallyObject and set the value ! to its default, as-initialized values, including deallocations. - this % label = "" + this % name = "" if (allocated(this % filters)) then do i = 1, size(this % filters) diff --git a/src/utils/openmc/material.py b/src/utils/openmc/material.py index c7292416e6..44277471b4 100644 --- a/src/utils/openmc/material.py +++ b/src/utils/openmc/material.py @@ -383,6 +383,9 @@ class Material(object): element = ET.Element("material") element.set("id", str(self._id)) + if len(self._name) > 0: + element.set("name", str(self._name)) + # Create density XML subelement subelement = ET.SubElement(element, "density") if self._density_units is not 'sum': @@ -520,10 +523,6 @@ class MaterialsFile(object): for material in self._materials: xml_element = material.get_material_xml() - - if len(material._name) > 0: - self._materials_file.append(ET.Comment(material._name)) - self._materials_file.append(xml_element) diff --git a/src/utils/particle_restart.py b/src/utils/openmc/particle_restart.py similarity index 100% rename from src/utils/particle_restart.py rename to src/utils/openmc/particle_restart.py diff --git a/src/utils/openmc/statepoint.py b/src/utils/openmc/statepoint.py index 73ec85fc2c..b641365de9 100644 --- a/src/utils/openmc/statepoint.py +++ b/src/utils/openmc/statepoint.py @@ -1,5 +1,6 @@ import copy import struct +import sys import numpy as np import scipy.stats @@ -7,6 +8,9 @@ import scipy.stats import openmc from openmc.constants import * +if sys.version > '3': + long = int + class SourceSite(object): @@ -111,6 +115,17 @@ class StatePoint(object): def tallies(self): return self._tallies + @property + def tallies_present(self): + return self._tallies_present + + @property + def global_tallies(self): + return self._global_tallies + + @property + def n_realizations(self): + return self._n_realizations def _read_metadata(self): @@ -298,18 +313,6 @@ class StatePoint(object): # Iterate over all Tallies for tally_key in self._tally_keys: - # Read user-specified Tally label (if specified) - label_size = self._get_int( - path='{0}{1}/label_size'.format(base, tally_key))[0] - - if label_size > 0: - label = self._get_string( - label_size, path='{0}{1}/label'.format(base, tally_key)) - - # Remove leading and trailing characters from string label - label = label.lstrip('[\'') - label = label.rstrip('\']') - # Read integer Tally estimator type code (analog or tracklength) estimator_type = self._get_int( path='{0}{1}/estimator'.format(base, tally_key))[0] @@ -319,7 +322,7 @@ class StatePoint(object): path='{0}{1}/n_realizations'.format(base, tally_key))[0] # Create Tally object and assign basic properties - tally = openmc.Tally(tally_key, label) + tally = openmc.Tally(tally_key) tally.estimator = ESTIMATOR_TYPES[estimator_type] tally.num_realizations = n_realizations @@ -367,7 +370,8 @@ class StatePoint(object): filter.num_bins = n_bins if FILTER_TYPES[filter_type] == 'mesh': - filter.mesh = self._meshes[bins] + key = self._mesh_keys[self._mesh_ids.index(bins)] + filter.mesh = self._meshes[key] # Add Filter to the Tally tally.add_filter(filter) @@ -570,7 +574,7 @@ class StatePoint(object): def get_tally(self, score, filters, nuclides, - label='', estimator='tracklength'): + name='', estimator='tracklength'): """Finds and returns a Tally object with certain properties. Parameters @@ -584,8 +588,8 @@ class StatePoint(object): nuclides : list A list of Nuclide objects - label : str - The label specified for the Tally (default is '') + name : str + The name specified for the Tally (default is '') estimator: str The type of estimator ('tracklength' (default) or 'analog') @@ -597,8 +601,8 @@ class StatePoint(object): # Iterate over all tallies to find the appropriate one for tally_id, test_tally in self._tallies.items(): - # Determine if the queried Tally label is the same as this Tally - if not label == test_tally._label: + # Determine if the queried Tally name is the same as this Tally + if not name == test_tally._name: continue # Determine if the queried Tally estimator is the same as this Tally @@ -609,7 +613,7 @@ class StatePoint(object): if not score in test_tally._scores: continue - # Determine if queried Tally filters is the same length as this Tally + # Determine if queried Tally filters is same length as this Tally if len(filters) != len(test_tally._filters): continue @@ -643,7 +647,7 @@ class StatePoint(object): return tally - def get_tally_id(self, score, filters, label='', estimator='tracklength'): + def get_tally_id(self, score, filters, name='', estimator='tracklength'): """Retrieve the Tally ID for a given list of filters and score(s). Parameters @@ -654,14 +658,14 @@ class StatePoint(object): filters : list A list of Filter objects - label : str - The label specified for the Tally (default is '') + name : str + The name specified for the Tally (default is '') estimator: str The type of estimator ('tracklength' (default) or 'analog') """ - tally = self.get_tally(score, filters, label, estimator) + tally = self.get_tally(score, filters, name, estimator) return tally._id @@ -674,6 +678,9 @@ class StatePoint(object): for tally_id, tally in self._tallies.items(): + # Get the Tally name from the summary file + tally.name = summary.tallies[tally_id].name + nuclide_zaids = copy.deepcopy(tally._nuclides) for nuclide_zaid in nuclide_zaids: diff --git a/src/utils/openmc/summary.py b/src/utils/openmc/summary.py index 0e74a28720..8aa25ad90e 100644 --- a/src/utils/openmc/summary.py +++ b/src/utils/openmc/summary.py @@ -26,6 +26,7 @@ class Summary(object): self._read_metadata() self._read_geometry() + self._read_tallies() def _read_metadata(self): @@ -105,6 +106,7 @@ class Summary(object): material_id = int(key.lstrip('material ')) index = self._f['materials'][key]['index'][0] + name = self._f['materials'][key]['name'][0] density = self._f['materials'][key]['atom_density'][0] nuc_densities = self._f['materials'][key]['nuclide_densities'][...] nuclides = self._f['materials'][key]['nuclides'][...] @@ -123,7 +125,7 @@ class Summary(object): sab_xs.append(sab_table.split('.')[1]) # Create the Material - material = openmc.Material(material_id=material_id) + material = openmc.Material(material_id=material_id, name=name) # Set the Material's density to g/cm3 - this is what is used in OpenMC material.set_density(density=density, units='g/cm3') @@ -164,6 +166,7 @@ class Summary(object): surface_id = int(key.lstrip('surface ')) index = self._f['geometry/surfaces'][key]['index'][0] + name = self._f['geometry/surfaces'][key]['name'][0] surf_type = self._f['geometry/surfaces'][key]['type'][...][0] bc = self._f['geometry/surfaces'][key]['boundary_condition'][...][0] coeffs = self._f['geometry/surfaces'][key]['coefficients'][...] @@ -172,47 +175,47 @@ class Summary(object): if surf_type == 'X Plane': x0 = coeffs[0] - surface = openmc.XPlane(surface_id, bc, x0) + surface = openmc.XPlane(surface_id, bc, x0, name) elif surf_type == 'Y Plane': y0 = coeffs[0] - surface = openmc.YPlane(surface_id, bc, y0) + surface = openmc.YPlane(surface_id, bc, y0, name) elif surf_type == 'Z Plane': z0 = coeffs[0] - surface = openmc.ZPlane(surface_id, bc, z0) + surface = openmc.ZPlane(surface_id, bc, z0, name) elif surf_type == 'Plane': A = coeffs[0] B = coeffs[1] C = coeffs[2] D = coeffs[3] - surface = openmc.Plane(surface_id, bc, A, B, C, D) + surface = openmc.Plane(surface_id, bc, A, B, C, D, name) elif surf_type == 'X Cylinder': y0 = coeffs[0] z0 = coeffs[1] R = coeffs[2] - surface = openmc.XCylinder(surface_id, bc, y0, z0, R) + surface = openmc.XCylinder(surface_id, bc, y0, z0, R, name) elif surf_type == 'Y Cylinder': x0 = coeffs[0] z0 = coeffs[1] R = coeffs[2] - surface = openmc.YCylinder(surface_id, bc, x0, z0, R) + surface = openmc.YCylinder(surface_id, bc, x0, z0, R, name) elif surf_type == 'Z Cylinder': x0 = coeffs[0] y0 = coeffs[1] R = coeffs[2] - surface = openmc.ZCylinder(surface_id, bc, x0, y0, R) + surface = openmc.ZCylinder(surface_id, bc, x0, y0, R, name) elif surf_type == 'Sphere': x0 = coeffs[0] y0 = coeffs[1] z0 = coeffs[2] R = coeffs[3] - surface = openmc.Sphere(surface_id, bc, x0, y0, z0, R) + surface = openmc.Sphere(surface_id, bc, x0, y0, z0, R, name) elif surf_type in ['X Cone', 'Y Cone', 'Z Cone']: x0 = coeffs[0] @@ -221,11 +224,11 @@ class Summary(object): R2 = coeffs[3] if surf_type == 'X Cone': - surface = openmc.XCone(surface_id, bc, x0, y0, z0, R2) + surface = openmc.XCone(surface_id, bc, x0, y0, z0, R2, name) if surf_type == 'Y Cone': - surface = openmc.YCone(surface_id, bc, x0, y0, z0, R2) + surface = openmc.YCone(surface_id, bc, x0, y0, z0, R2, name) if surf_type == 'Z Cone': - surface = openmc.ZCone(surface_id, bc, x0, y0, z0, R2) + surface = openmc.ZCone(surface_id, bc, x0, y0, z0, R2, name) # Add Surface to global dictionary of all Surfaces self.surfaces[index] = surface @@ -255,6 +258,7 @@ class Summary(object): cell_id = int(key.lstrip('cell ')) index = self._f['geometry/cells'][key]['index'][0] + name = self._f['geometry/cells'][key]['name'][0] fill_type = self._f['geometry/cells'][key]['fill_type'][...][0] if fill_type == 'normal': @@ -270,7 +274,7 @@ class Summary(object): surfaces = [] # Create this Cell - cell = openmc.Cell(cell_id=cell_id) + cell = openmc.Cell(cell_id=cell_id, name=name) if fill_type == 'universe': translated = self._f['geometry/cells'][key]['translated'][0] @@ -348,6 +352,7 @@ class Summary(object): lattice_id = int(key.lstrip('lattice ')) index = self._f['geometry/lattices'][key]['index'][0] + name = self._f['geometry/lattices'][key]['name'][0] lattice_type = self._f['geometry/lattices'][key]['type'][...][0] if lattice_type == 'rectangular': @@ -363,11 +368,12 @@ class Summary(object): universe_ids = np.swapaxes(universe_ids, 1, 2) # Create the Lattice - lattice = openmc.RectLattice(lattice_id=lattice_id) + lattice = openmc.RectLattice(lattice_id=lattice_id, name=name) lattice.dimension = tuple(dimension) lattice.lower_left = lower_left lattice.pitch = pitch + # If the Universe specified outer the Lattice is not void (-22) if outer != -22: lattice.outer = self.universes[outer] @@ -393,8 +399,8 @@ class Summary(object): self.lattices[index] = lattice if lattice_type == 'hexagonal': - n_rings = self._f['geometry/latties'][key]['n_rings'][0] - n_axial = self._f['geometry/latties'][key]['n_axial'][0] + n_rings = self._f['geometry/lattices'][key]['n_rings'][0] + n_axial = self._f['geometry/lattices'][key]['n_axial'][0] center = self._f['geometry/lattices'][key]['center'][...] pitch = self._f['geometry/lattices'][key]['pitch'][...] outer = self._f['geometry/lattices'][key]['outer'][0] @@ -403,8 +409,8 @@ class Summary(object): self._f['geometry/lattices'][key]['universes'][...] # Create the Lattice - lattice = openmc.HexLattice(lattice_id=lattice_id) - lattice.num_rings(n_rings) + lattice = openmc.HexLattice(lattice_id=lattice_id, name=name) + lattice.num_rings = n_rings lattice.num_axial = n_axial lattice.center = center lattice.pitch = pitch @@ -459,6 +465,73 @@ class Summary(object): self.openmc_geometry.root_universe = root_universe + def _read_tallies(self): + + # Initialize dictionaries for the Tallies + # Keys - Tally IDs + # Values - Tally objects + self.tallies = {} + + # Read the number of tallies + self.n_tallies = self._f['tallies/n_tallies'][0] + + # OpenMC Tally keys + all_keys = self._f['tallies/'].keys() + tally_keys = [key for key in all_keys if 'tally' in key] + + base = 'tallies/tally ' + + # Iterate over all Tallies + for tally_key in tally_keys: + + tally_id = int(tally_key.strip('tally ')) + subbase = '{0}{1}'.format(base, tally_id) + + # Read Tally name metadata + name_size = self._f['{0}/name_size'.format(subbase)][0] + if (name_size > 0): + tally_name = self._f['{0}/name'.format(subbase)][0] + tally_name = tally_name.lstrip('[\'') + tally_name = tally_name.rstrip('\']') + else: + tally_name = '' + + # Create Tally object and assign basic properties + tally = openmc.Tally(tally_id, tally_name) + + # Read score metadata + score_bins = self._f['{0}/score_bins'.format(subbase)][...] + for score_bin in score_bins: + tally.add_score(openmc.SCORE_TYPES[score_bin]) + num_score_bins = self._f['{0}/n_score_bins'.format(subbase)][...] + tally.num_score_bins = num_score_bins + + # Read filter metadata + num_filters = self._f['{0}/n_filters'.format(subbase)][0] + + # Initialize all Filters + for j in range(1, num_filters+1): + + subsubbase = '{0}/filter {1}'.format(subbase, j) + + # Read filter type (e.g., "cell", "energy", etc.) + filter_type = self._f['{0}/type_name'.format(subsubbase)][0] + + # Read the filter bins + num_bins = self._f['{0}/n_bins'.format(subsubbase)][0] + bins = self._f['{0}/bins'.format(subsubbase)][...] + + # Create Filter object + filter = openmc.Filter(filter_type, bins) + filter.num_bins = num_bins + + # Add Filter to the Tally + tally.add_filter(filter) + + # Add Tally to the global dictionary of all Tallies + self.tallies[tally_id] = tally + + def make_opencg_geometry(self): if self.opencg_geometry is None: self.opencg_geometry = get_opencg_geometry(self.openmc_geometry) diff --git a/src/utils/openmc/surface.py b/src/utils/openmc/surface.py index ab7ad5a7df..0e958b4565 100644 --- a/src/utils/openmc/surface.py +++ b/src/utils/openmc/surface.py @@ -133,6 +133,10 @@ class Surface(object): element = ET.Element("surface") element.set("id", str(self._id)) + + if len(self._name) > 0: + element.set("name", str(self._name)) + element.set("type", self._type) element.set("boundary", self._boundary_type) diff --git a/src/utils/openmc/tallies.py b/src/utils/openmc/tallies.py index 810611b26c..1dc8e536e0 100644 --- a/src/utils/openmc/tallies.py +++ b/src/utils/openmc/tallies.py @@ -584,11 +584,11 @@ class Mesh(object): class Tally(object): - def __init__(self, tally_id=None, label=''): + def __init__(self, tally_id=None, name=''): # Initialize Tally class attributes self.id = tally_id - self.label = label + self.name = name self._filters = [] self._nuclides = [] self._scores = [] @@ -612,7 +612,7 @@ class Tally(object): clone = type(self).__new__(type(self)) clone._id = self._id - clone._label = self._label + clone._name = self._name clone._estimator = self._estimator clone._num_score_bins = self._num_score_bins clone._num_realizations = self._num_realizations @@ -678,7 +678,7 @@ class Tally(object): hashable.append(score) hashable.append(self._estimator) - hashable.append(self._label) + hashable.append(self._name) return hash(tuple(hashable)) @@ -699,8 +699,8 @@ class Tally(object): @property - def label(self): - return self._label + def name(self): + return self._name @property @@ -819,16 +819,16 @@ class Tally(object): self._id = tally_id - @label.setter - def label(self, label): + @name.setter + def name(self, name): - if not is_string(label): + if not is_string(name): msg = 'Unable to set name for Tally ID={0} with a non-string ' \ - 'value {1}'.format(self._id, label) + 'value {1}'.format(self._id, name) raise ValueError(msg) else: - self._label = label + self._name = name def add_filter(self, filter): @@ -945,7 +945,7 @@ class Tally(object): string = 'Tally\n' string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) - string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._label) + string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) string += '{0: <16}\n'.format('\tFilters') @@ -976,9 +976,9 @@ class Tally(object): # Tally ID element.set("id", str(self._id)) - # Optional Tally label - if self._label != '': - element.set("label", self._label) + # Optional Tally name + if self._name != '': + element.set("name", self._name) # Optional Tally filters for filter in self._filters: @@ -1124,10 +1124,10 @@ class Tally(object): # Convert (x,y,z) to a single bin -- this is similar to # subroutine mesh_indices_to_bin in openmc/src/mesh.F90. - value = ((filter_bins[i][0] - 1) * ny * nz + - (filter_bins[i][1] - 1) * nz + - (filter_bins[i][2] - 1)) - filter_index += value * test_filter._stride + val = ((filter_bins[i][0] - 1) * ny * nz + + (filter_bins[i][1] - 1) * nz + + (filter_bins[i][2] - 1)) + filter_index += val * test_filter._stride # Filter bins for distribcell are the "IDs" of each unique placement # of the Cell in the Geometry (integers starting at 0) @@ -1220,7 +1220,7 @@ class Tally(object): # Add basic Tally data to the HDF5 group tally_group.create_dataset('id', data=self._id) - tally_group.create_dataset('label', data=self._label) + tally_group.create_dataset('name', data=self._name) tally_group.create_dataset('estimator', data=self._estimator) tally_group.create_dataset('scores', data=np.array(self._scores)) @@ -1268,7 +1268,7 @@ class Tally(object): # Add basic Tally data to the nested dictionary tally_group['id'] = self._id - tally_group['label'] = self._label + tally_group['name'] = self._name tally_group['estimator'] = self._estimator tally_group['scores'] = np.array(self._scores) diff --git a/src/utils/openmc/universe.py b/src/utils/openmc/universe.py index ab663269da..ee05e5356d 100644 --- a/src/utils/openmc/universe.py +++ b/src/utils/openmc/universe.py @@ -345,6 +345,9 @@ class Cell(object): element = ET.Element("cell") element.set("id", str(self._id)) + if len(self._name) > 0: + element.set("name", str(self._name)) + if isinstance(self._fill, openmc.Material): element.set("material", str(self._fill._id)) @@ -376,10 +379,6 @@ class Cell(object): # Create the XML subelement for this Surface surface = self._surfaces[surface_id][0] surface_subelement = surface.create_xml_subelement() - - if len(surface._name) > 0: - xml_element.append(ET.Comment(surface._name)) - xml_element.append(surface_subelement) # Append the halfspace and Surface ID @@ -614,10 +613,6 @@ class Universe(object): # Append the Universe ID to the subelement and add to Element cell_subelement.set("universe", str(self._id)) - - if len(cell._name) > 0: - xml_element.append(ET.Comment(cell._name)) - xml_element.append(cell_subelement) @@ -995,6 +990,9 @@ class RectLattice(Lattice): lattice_subelement = ET.Element("lattice") lattice_subelement.set("id", str(self._id)) + if len(self._name) > 0: + lattice_subelement.set("name", str(self._name)) + # Export the Lattice cell pitch if len(self._pitch) == 3: pitch = ET.SubElement(lattice_subelement, "pitch") @@ -1078,9 +1076,6 @@ class RectLattice(Lattice): universes = ET.SubElement(lattice_subelement, "universes") universes.text = universe_ids - if len(self._name) > 0: - xml_element.append(ET.Comment(self._name)) - # Append the XML subelement for this Lattice to the XML element xml_element.append(lattice_subelement) @@ -1324,6 +1319,9 @@ class HexLattice(Lattice): lattice_subelement = ET.Element("hex_lattice") lattice_subelement.set("id", str(self._id)) + if len(self._name) > 0: + lattice_subelement.set("name", str(self._name)) + # Export the Lattice cell pitch if len(self._pitch) == 2: pitch = ET.SubElement(lattice_subelement, "pitch") @@ -1394,9 +1392,6 @@ class HexLattice(Lattice): universes = ET.SubElement(lattice_subelement, "universes") universes.text = '\n' + universe_ids - if len(self._name) > 0: - xml_element.append(ET.Comment(self._name)) - # Append the XML subelement for this Lattice to the XML element xml_element.append(lattice_subelement) diff --git a/src/utils/plot_mesh_tally.py b/src/utils/plot_mesh_tally.py index b78305ebb6..5c46d927c0 100755 --- a/src/utils/plot_mesh_tally.py +++ b/src/utils/plot_mesh_tally.py @@ -11,7 +11,7 @@ from matplotlib.figure import Figure import matplotlib.pyplot as plt import numpy as np -from statepoint import * +from openmc.statepoint import StatePoint if sys.version_info[0] < 3: import Tkinter as tk @@ -130,14 +130,13 @@ class MeshPlotter(tk.Frame): selectedTally = self.datafile.tallies[tally_id] # Get mesh for selected tally - self.mesh = self.datafile.meshes[ - selectedTally.filters['mesh'].bins[0] - 1] + self.mesh = selectedTally.filters_by_name['mesh'].mesh # Get mesh dimensions self.nx, self.ny, self.nz = self.mesh.dimension # Repopulate comboboxes baesd on current basis selection - text = self.basisBox['values'][self.basisBox.current()] + text = self.basisBox.get() if text == 'xy': self.axialBox['values'] = [str(i+1) for i in range(self.nz)] elif text == 'yz': @@ -164,7 +163,8 @@ class MeshPlotter(tk.Frame): # create a label/combobox for each filter in selected tally count = 0 - for filterType in selectedTally.filters: + for f in selectedTally.filters: + filterType = f.type if filterType == 'mesh': continue count += 1 @@ -176,7 +176,6 @@ class MeshPlotter(tk.Frame): self.filterBoxes[filterType] = combobox # Set combobox items - f = selectedTally.filters[filterType] if filterType in ['energyin', 'energyout']: combobox['values'] = ['{0} to {1}'.format(*f.bins[i:i+2]) for i in range(f.length)] @@ -198,7 +197,7 @@ class MeshPlotter(tk.Frame): def redraw(self, event=None): basis = self.basisBox.current() + 1 axial_level = self.axialBox.current() + 1 - is_mean = self.meanBox.current() + mbvalue = self.meanBox.get() # Get selected tally tally_id = self.meshTallies[self.tallyBox.current()] @@ -206,71 +205,45 @@ class MeshPlotter(tk.Frame): # Create spec_list spec_list = [] - for f in selectedTally.filters.values(): + for f in selectedTally.filters: if f.type == 'mesh': + mesh_filter = f continue index = self.filterBoxes[f.type].current() - spec_list.append((f.type, index)) + spec_list.append((f, index)) - # Take is_mean and convert it to an index of the score - score_loc = is_mean - if score_loc > 1: - score_loc = 1 - - text = self.basisBox['values'][self.basisBox.current()] + text = self.basisBox.get() if text == 'xy': - matrix = np.zeros((self.nx, self.ny)) - for i in range(self.nx): - for j in range(self.ny): - matrix[i, j] = self.datafile.get_value(tally_id, - spec_list + [('mesh', (i + 1, j + 1, axial_level))], - self.scoreBox.current())[score_loc] - # Calculate relative uncertainty from absolute, if requested - if is_mean == 2: - # Take care to handle zero means when normalizing - mean_val = self.datafile.get_value(tally_id, - spec_list + [('mesh', (i + 1, j + 1, axial_level))], - self.scoreBox.current())[0] - if mean_val > 0.0: - matrix[i, j] = matrix[i, j] / mean_val - else: - matrix[i, j] = 0.0 - + dims = (self.nx, self.ny) elif text == 'yz': - matrix = np.zeros((self.ny, self.nz)) - for i in range(self.ny): - for j in range(self.nz): - matrix[i, j] = self.datafile.get_value(tally_id, - spec_list + [('mesh', (axial_level, i + 1, j + 1))], - self.scoreBox.current())[score_loc] - # Calculate relative uncertainty from absolute, if requested - if is_mean == 2: - # Take care to handle zero means when normalizing - mean_val = self.datafile.get_value(tally_id, - spec_list + [('mesh', (axial_level, i + 1, j + 1))], - self.scoreBox.current())[0] - if mean_val > 0.0: - matrix[i, j] = matrix[i, j] / mean_val - else: - matrix[i, j] = 0.0 - + dims = (self.ny, self.nz) else: - matrix = np.zeros((self.nx, self.nz)) - for i in range(self.nx): - for j in range(self.nz): - matrix[i, j] = self.datafile.get_value(tally_id, - spec_list + [('mesh', (i + 1, axial_level, j + 1))], - self.scoreBox.current())[score_loc] - # Calculate relative uncertainty from absolute, if requested - if is_mean == 2: - # Take care to handle zero means when normalizing - mean_val = self.datafile.get_value(tally_id, - spec_list + [('mesh', (i + 1, axial_level, j + 1))], - self.scoreBox.current())[0] - if mean_val > 0.0: - matrix[i, j] = matrix[i, j] / mean_val - else: - matrix[i, j] = 0.0 + dims = (self.nx, self.nz) + matrix = np.zeros(dims) + for i in range(dims[0]): + for j in range(dims[1]): + if text == 'xy': + meshtuple = (i + 1, j + 1, axial_level) + elif text == 'yz': + meshtuple = (axial_level, i + 1, j + 1) + else: + meshtuple = (i + 1, axial_level, j + 1) + filters, filter_bins = zip(*spec_list + [ + (mesh_filter, meshtuple)]) + mean = selectedTally.get_value( + self.scoreBox.get(), filters, filter_bins) + stdev = selectedTally.get_value( + self.scoreBox.get(), filters, filter_bins, + value='std_dev') + if mbvalue == 'Mean': + matrix[i, j] = mean + elif mbvalue == 'Absolute uncertainty': + matrix[i, j] = stdev + else: + if mean > 0.: + matrix[i, j] = stdev/mean + else: + matrix[i, j] = 0. # Clear the figure self.fig.clear() @@ -292,13 +265,14 @@ class MeshPlotter(tk.Frame): # Create StatePoint object and read in data self.datafile = StatePoint(filename) self.datafile.read_results() - self.datafile.generate_stdev() + self.datafile.compute_stdev() # Find which tallies are mesh tallies self.meshTallies = [] - for itally, tally in enumerate(self.datafile.tallies): - if 'mesh' in tally.filters: + for itally, tally in self.datafile.tallies.iteritems(): + if any([f.type == 'mesh' for f in tally.filters]): self.meshTallies.append(itally) + tally.filters_by_name = {f.type: f for f in tally.filters} if not self.meshTallies: messagebox.showerror("Invalid StatePoint File", diff --git a/src/utils/statepoint_cmp.py b/src/utils/statepoint_cmp.py deleted file mode 100755 index 7ad284370c..0000000000 --- a/src/utils/statepoint_cmp.py +++ /dev/null @@ -1,86 +0,0 @@ -#!/usr/bin/env python - -import sys - -from numpy.testing import assert_allclose, assert_equal - -from openmc.statepoint import StatePoint - -if len(sys.argv) > 2: - path1 = sys.argv[1] - path2 = sys.argv[2] -else: - raise - -# Create StatePoint objects -sp1 = StatePoint(path1) -sp2 = StatePoint(path2) - -# Read tally results -sp1.read_results() -sp2.read_results() - -# Compare header information -assert sp1.revision == sp2.revision -assert sp1.version == sp2.version -assert sp1.seed == sp2.seed -assert sp1.run_mode == sp2.run_mode -assert sp1.n_particles == sp2.n_particles -assert sp1.n_batches == sp2.n_batches - -if sp1.run_mode == 2: - # Compare criticality information - assert sp1.n_inactive == sp2.n_inactive - assert sp1.gen_per_batch == sp2.gen_per_batch - assert sp1.current_batch == sp2.current_batch - - # Compare keff results - assert_allclose(sp1.k_batch, sp2.k_batch) - - # Compare entropy results - assert_allclose(sp1.entropy, sp2.entropy) - -# Compare global tallies -assert_allclose(sp1.global_tallies, sp2.global_tallies) - -# Compare meshes -assert len(sp1.meshes) == len(sp2.meshes) -for m1, m2 in zip(sp1.meshes, sp2.meshes): - assert m1.type == m2.type - assert_equal(m1.dimension, m2.dimension) - assert_equal(m1.lower_left, m2.lower_left) - assert_equal(m1.upper_right, m2.upper_right) - assert_equal(m1.width, m2.width) - -# Compare tallies -assert len(sp1.tallies) == len(sp2.tallies) -for t1, t2 in zip(sp1.tallies, sp2.tallies): - # Compare size of tallies - assert t1.total_score_bins == t2.total_score_bins - assert t1.total_filter_bins == t2.total_filter_bins - - # Compare filters - assert len(t1.filters) == len(t2.filters) - for f1, f2 in zip(t1.filters.values(), t2.filters.values()): - assert f1.type == f2.type - assert f1.length == f2.length - assert_equal(f1.bins, f2.bins) - - # Compare nuclide and score bins - assert_equal(t1.nuclides, t2.nuclides) - assert_equal(t1.scores, t2.scores) - - # Compare tally results - assert_allclose(t1.results, t2.results) - -# If criticality, compare source sites -if sp1.run_mode == 2: - sp1.read_source() - sp2.read_source() - - assert len(sp1.source) == len(sp2.source) - for s1, s2 in zip(sp1.source, sp2.source): - assert s1.weight == s2.weight - assert_allclose(s1.xyz, s2.xyz) - assert_allclose(s1.uvw, s2.uvw) - assert s1.E == s2.E diff --git a/src/utils/statepoint_histogram.py b/src/utils/statepoint_histogram.py index ed7b2682b9..26e8bdae60 100755 --- a/src/utils/statepoint_histogram.py +++ b/src/utils/statepoint_histogram.py @@ -13,60 +13,31 @@ from openmc.statepoint import StatePoint # Get filename filename = argv[1] -# Determine axial level -axial_level = int(argv[2]) if len(argv) > 2 else 1 -score = int(argv[3]) if len(argv) > 3 else 1 - # Create StatePoint object sp = StatePoint(filename) - -# Read number of realizations for global tallies -sp.n_realizations = sp._get_int()[0] - -# Read global tallies -n_global_tallies = sp._get_int()[0] -sp.global_tallies = np.array(sp._get_double(2*n_global_tallies)) -sp.global_tallies.shape = (n_global_tallies, 2) - -# Flag indicating if tallies are present -tallies_present = sp._get_int()[0] +sp.read_results() +sp.compute_ci() # Check if tallies are present -if not tallies_present: +if not sp.tallies_present: raise Exception("No tally data in state point!") # Loop over all tallies -for i, t in enumerate(sp.tallies): - # Calculate t-value for 95% two-sided CI - n = t.n_realizations - t_value = scipy.stats.t.ppf(0.975, n - 1) +for i, t in sp.tallies.items(): + # Determine relative error and fraction of bins with less than 1% half-width + # of CI + n_bins = t.mean.size + relative_error = t.std_dev[t.mean > 0.] / t.mean[t.mean > 0.] + fraction = float(sum(relative_error < 0.01))/n_bins - # Create lists for tallies - mean = [] - uncertainties = [] - nonzero = [] - - n_bins = t.total_score_bins * t.total_filter_bins - - # Loop over filter/score bins - for j in range(n_bins): - # Read sum and sum-squared - s, s2 = sp._get_double(2) - s /= n - if s != 0.0: - relative_error = t_value*sqrt((s2/n - s*s)/(n-1))/s - uncertainties.append(relative_error) - - # Display number of tallies with less than 1% CIs - fraction = (sum([1.0 if re < 0.01 else 0.0 for re in uncertainties]) - / n_bins) - print("Tally " + str(i+1)) + # Display results + print("Tally " + str(i)) print(" Fraction under 1% = {0}".format(fraction)) - print(" Min relative error = {0}".format(min(uncertainties))) - print(" Max relative error = {0}".format(max(uncertainties))) + print(" Min relative error = {0}".format(min(relative_error))) + print(" Max relative error = {0}".format(max(relative_error))) print(" Non-scoring bins = {0}".format( - 1.0 - float(len(uncertainties))/n_bins)) + 1.0 - float(relative_error.size)/n_bins)) # Plot histogram - plt.hist(uncertainties, 100) + plt.hist(relative_error, 100) plt.show() diff --git a/src/utils/statepoint_meshplot.py b/src/utils/statepoint_meshplot.py deleted file mode 100755 index ba09ef5e68..0000000000 --- a/src/utils/statepoint_meshplot.py +++ /dev/null @@ -1,95 +0,0 @@ -#!/usr/bin/env python - -from __future__ import print_function, division - -from sys import argv -from math import sqrt - -import numpy as np -import scipy.stats -import matplotlib.pyplot as plt - -from openmc.statepoint import StatePoint - -# Get filename -filename = argv[1] - -# Determine axial level -axial_level = int(argv[2]) if len(argv) > 2 else 1 -score = int(argv[3]) if len(argv) > 3 else 1 - -# Create StatePoint object -sp = StatePoint(filename) - -# Read number of realizations for global tallies -sp.n_realizations = sp._get_int()[0] - -# Read global tallies -n_global_tallies = sp._get_int()[0] -sp.global_tallies = np.array(sp._get_double(2*n_global_tallies)) -sp.global_tallies.shape = (n_global_tallies, 2) - -# Flag indicating if tallies are present -tallies_present = sp._get_int()[0] - -# Check if tallies are present -if not tallies_present: - raise Exception("No tally data in state point!") - -# Loop over all tallies -print("Reading data...") -for t in sp.tallies: - # Calculate t-value for 95% two-sided CI - n = t.n_realizations - t_value = scipy.stats.t.ppf(0.975, n - 1) - - n_bins = t.total_score_bins * t.total_filter_bins - - # Check for mesh - if 'mesh' not in t.filters: - continue - - # Get Mesh object and determine size - m = sp.meshes[t.filters['mesh'].bins[0] - 1] - if len(m.dimension) == 2: - nx, ny = m.dimension - nz = 1 - else: - nx, ny, nz = m.dimension - - # Calculate number of score bins - ns = t.total_score_bins * t.total_filter_bins // (nx*ny*nz) - assert n_bins == nx*ny*nz*ns - - # Create lists for tallies - mean = np.zeros((nx, ny)) - error = np.zeros((nx, ny)) - criteria = np.zeros((nx, ny)) - - # Determine starting position for data - start = sp._f.tell() + (axial_level-1)*ns*16 + (score - 1)*16 - - for x in range(nx): - for y in range(ny): - # Seek to position of data - sp._f.seek(start + x*ny*nz*ns*16 + y*nz*ns*16) - - # Read sum and sum-squared - s, s2 = sp._get_double(2) - s /= n - mean[x, y] = s - if s != 0.0: - error[x, y] = t_value*sqrt((s2/n - s*s)/(n-1))/s - criteria[x, y] = 1.0 if error[x, y] < 0.05 else 0.0 - - # Make figure - print("Making colorplot...") - plt.imshow(mean, interpolation="nearest") - plt.colorbar() - plt.xlim((0, nx)) - plt.ylim((0, ny)) - plt.xticks(np.linspace(0, nx, 5), - np.linspace(m.lower_left[0], m.upper_right[0], 5)) - plt.yticks(np.linspace(0, ny, 5), - np.linspace(m.lower_left[1], m.upper_right[1], 5)) - plt.show() diff --git a/tests/test_cmfd_feed/results.py b/tests/test_cmfd_feed/results.py index b55ef81311..125d13a8ea 100644 --- a/tests/test_cmfd_feed/results.py +++ b/tests/test_cmfd_feed/results.py @@ -6,6 +6,7 @@ import numpy as np sys.path.insert(0, '../../src/utils') # import statepoint +import openmc from openmc.statepoint import StatePoint # read in statepoint file @@ -17,37 +18,38 @@ else: sp.read_results() # extract tally results and convert to vector -tally1 = sp._tallies[1] -results1 = np.zeros((tally1._sum.size + tally1._sum.size, )) -results1[0::2] = tally1._sum.ravel() -results1[1::2] = tally1._sum_sq.ravel() +tally1 = sp.get_tally('flux', [openmc.Filter('mesh', [1])], \ + [-1], estimator='tracklength') +tally2 = sp.get_tally('flux', [openmc.Filter('mesh', [2])], \ + [-1], estimator='analog') +tally3 = sp.get_tally('nu-fission', [openmc.Filter('mesh', [2])], \ + [-1], estimator='analog') +tally4 = sp.get_tally('current', [openmc.Filter('mesh', [2]), \ + openmc.Filter('surface', [1,2,3,4,5,6])], \ + [-1], estimator='analog') -for tally_id in sp._tallies: - if 'CMFD flux, total, scatter-1' in sp._tallies[tally_id]._label: - tally2 = sp._tallies[tally_id] - elif 'CMFD neutron production' in sp._tallies[tally_id]._label: - tally3 = sp._tallies[tally_id] - elif 'CMFD surface currents' in sp._tallies[tally_id]._label: - tally4 = sp._tallies[tally_id] +results1 = np.zeros((tally1.sum.size + tally1.sum.size, )) +results1[0::2] = tally1.sum.ravel() +results1[1::2] = tally1.sum_sq.ravel() -results2 = np.zeros((tally2._sum.size + tally2._sum.size, )) -results2[0::2] = tally2._sum.ravel() -results2[1::2] = tally2._sum_sq.ravel() +results2 = np.zeros((tally2.sum.size + tally2.sum.size, )) +results2[0::2] = tally2.sum.ravel() +results2[1::2] = tally2.sum_sq.ravel() -results3 = np.zeros((tally3._sum.size + tally3._sum.size, )) -results3[0::2] = tally3._sum.ravel() -results3[1::2] = tally3._sum_sq.ravel() +results3 = np.zeros((tally3.sum.size + tally3.sum.size, )) +results3[0::2] = tally3.sum.ravel() +results3[1::2] = tally3.sum_sq.ravel() -results4 = np.zeros((tally4._sum.size + tally4._sum.size, )) -results4[0::2] = tally4._sum.ravel() -results4[1::2] = tally4._sum_sq.ravel() +results4 = np.zeros((tally4.sum.size + tally4.sum.size, )) +results4[0::2] = tally4.sum.ravel() +results4[1::2] = tally4.sum_sq.ravel() # set up output string outstr = '' # write out k-combined outstr += 'k-combined:\n' -outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1]) +outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1]) # write out tally results outstr += 'tally 1:\n' diff --git a/tests/test_cmfd_jfnk/results.py b/tests/test_cmfd_jfnk/results.py index b55ef81311..e104fbd9ae 100644 --- a/tests/test_cmfd_jfnk/results.py +++ b/tests/test_cmfd_jfnk/results.py @@ -6,6 +6,7 @@ import numpy as np sys.path.insert(0, '../../src/utils') # import statepoint +import openmc from openmc.statepoint import StatePoint # read in statepoint file @@ -17,30 +18,31 @@ else: sp.read_results() # extract tally results and convert to vector -tally1 = sp._tallies[1] -results1 = np.zeros((tally1._sum.size + tally1._sum.size, )) -results1[0::2] = tally1._sum.ravel() -results1[1::2] = tally1._sum_sq.ravel() +tally1 = sp.get_tally('flux', [openmc.Filter('mesh', [1])], \ + [-1], estimator='tracklength') +tally2 = sp.get_tally('flux', [openmc.Filter('mesh', [2])], \ + [-1], estimator='analog') +tally3 = sp.get_tally('nu-fission', [openmc.Filter('mesh', [2])], \ + [-1], estimator='analog') +tally4 = sp.get_tally('current', [openmc.Filter('mesh', [2]), \ + openmc.Filter('surface', [1,2,3,4,5,6])], \ + [-1], estimator='analog') -for tally_id in sp._tallies: - if 'CMFD flux, total, scatter-1' in sp._tallies[tally_id]._label: - tally2 = sp._tallies[tally_id] - elif 'CMFD neutron production' in sp._tallies[tally_id]._label: - tally3 = sp._tallies[tally_id] - elif 'CMFD surface currents' in sp._tallies[tally_id]._label: - tally4 = sp._tallies[tally_id] +results1 = np.zeros((tally1.sum.size + tally1.sum.size, )) +results1[0::2] = tally1.sum.ravel() +results1[1::2] = tally1.sum_sq.ravel() -results2 = np.zeros((tally2._sum.size + tally2._sum.size, )) -results2[0::2] = tally2._sum.ravel() -results2[1::2] = tally2._sum_sq.ravel() +results2 = np.zeros((tally2.sum.size + tally2.sum.size, )) +results2[0::2] = tally2.sum.ravel() +results2[1::2] = tally2.sum_sq.ravel() -results3 = np.zeros((tally3._sum.size + tally3._sum.size, )) -results3[0::2] = tally3._sum.ravel() -results3[1::2] = tally3._sum_sq.ravel() +results3 = np.zeros((tally3.sum.size + tally3.sum.size, )) +results3[0::2] = tally3.sum.ravel() +results3[1::2] = tally3.sum_sq.ravel() -results4 = np.zeros((tally4._sum.size + tally4._sum.size, )) -results4[0::2] = tally4._sum.ravel() -results4[1::2] = tally4._sum_sq.ravel() +results4 = np.zeros((tally4.sum.size + tally4.sum.size, )) +results4[0::2] = tally4.sum.ravel() +results4[1::2] = tally4.sum_sq.ravel() # set up output string outstr = '' diff --git a/tests/test_cmfd_nofeed/results.py b/tests/test_cmfd_nofeed/results.py index b55ef81311..e104fbd9ae 100644 --- a/tests/test_cmfd_nofeed/results.py +++ b/tests/test_cmfd_nofeed/results.py @@ -6,6 +6,7 @@ import numpy as np sys.path.insert(0, '../../src/utils') # import statepoint +import openmc from openmc.statepoint import StatePoint # read in statepoint file @@ -17,30 +18,31 @@ else: sp.read_results() # extract tally results and convert to vector -tally1 = sp._tallies[1] -results1 = np.zeros((tally1._sum.size + tally1._sum.size, )) -results1[0::2] = tally1._sum.ravel() -results1[1::2] = tally1._sum_sq.ravel() +tally1 = sp.get_tally('flux', [openmc.Filter('mesh', [1])], \ + [-1], estimator='tracklength') +tally2 = sp.get_tally('flux', [openmc.Filter('mesh', [2])], \ + [-1], estimator='analog') +tally3 = sp.get_tally('nu-fission', [openmc.Filter('mesh', [2])], \ + [-1], estimator='analog') +tally4 = sp.get_tally('current', [openmc.Filter('mesh', [2]), \ + openmc.Filter('surface', [1,2,3,4,5,6])], \ + [-1], estimator='analog') -for tally_id in sp._tallies: - if 'CMFD flux, total, scatter-1' in sp._tallies[tally_id]._label: - tally2 = sp._tallies[tally_id] - elif 'CMFD neutron production' in sp._tallies[tally_id]._label: - tally3 = sp._tallies[tally_id] - elif 'CMFD surface currents' in sp._tallies[tally_id]._label: - tally4 = sp._tallies[tally_id] +results1 = np.zeros((tally1.sum.size + tally1.sum.size, )) +results1[0::2] = tally1.sum.ravel() +results1[1::2] = tally1.sum_sq.ravel() -results2 = np.zeros((tally2._sum.size + tally2._sum.size, )) -results2[0::2] = tally2._sum.ravel() -results2[1::2] = tally2._sum_sq.ravel() +results2 = np.zeros((tally2.sum.size + tally2.sum.size, )) +results2[0::2] = tally2.sum.ravel() +results2[1::2] = tally2.sum_sq.ravel() -results3 = np.zeros((tally3._sum.size + tally3._sum.size, )) -results3[0::2] = tally3._sum.ravel() -results3[1::2] = tally3._sum_sq.ravel() +results3 = np.zeros((tally3.sum.size + tally3.sum.size, )) +results3[0::2] = tally3.sum.ravel() +results3[1::2] = tally3.sum_sq.ravel() -results4 = np.zeros((tally4._sum.size + tally4._sum.size, )) -results4[0::2] = tally4._sum.ravel() -results4[1::2] = tally4._sum_sq.ravel() +results4 = np.zeros((tally4.sum.size + tally4.sum.size, )) +results4[0::2] = tally4.sum.ravel() +results4[1::2] = tally4.sum_sq.ravel() # set up output string outstr = '' diff --git a/tests/test_particle_restart_eigval/results.py b/tests/test_particle_restart_eigval/results.py index ae446b855d..2ff1ef2b89 100644 --- a/tests/test_particle_restart_eigval/results.py +++ b/tests/test_particle_restart_eigval/results.py @@ -5,7 +5,7 @@ import sys sys.path.insert(0, '../../src/utils') # import particle restart -import particle_restart as pr +import openmc.particle_restart as pr # read in particle restart file if len(sys.argv) > 1: diff --git a/tests/test_particle_restart_fixed/results.py b/tests/test_particle_restart_fixed/results.py index a5bac1da9b..6c9a40d5a8 100644 --- a/tests/test_particle_restart_fixed/results.py +++ b/tests/test_particle_restart_fixed/results.py @@ -4,8 +4,8 @@ import sys sys.path.insert(0, '../../src/utils') -# import particle restart -import particle_restart as pr +# import particle restart +import openmc.particle_restart as pr # read in particle restart file if len(sys.argv) > 1: @@ -15,8 +15,8 @@ else: # set up output string outstr = '' - -# write out properties + +# write out properties outstr += 'current batch:\n' outstr += "{0:12.6E}\n".format(p.current_batch) outstr += 'current gen:\n'