diff --git a/docs/img/3dcore.png b/docs/img/3dcore.png
new file mode 100644
index 0000000000..69fed4a079
Binary files /dev/null and b/docs/img/3dcore.png differ
diff --git a/docs/img/3dgeomplot.png b/docs/img/3dgeomplot.png
new file mode 100644
index 0000000000..5c16d8d831
Binary files /dev/null and b/docs/img/3dgeomplot.png differ
diff --git a/docs/img/atr.png b/docs/img/atr.png
new file mode 100644
index 0000000000..e1f440584b
Binary files /dev/null and b/docs/img/atr.png differ
diff --git a/docs/img/fluxplot.png b/docs/img/fluxplot.png
new file mode 100644
index 0000000000..9c13ff33f3
Binary files /dev/null and b/docs/img/fluxplot.png differ
diff --git a/docs/img/plotmeshtally.png b/docs/img/plotmeshtally.png
new file mode 100644
index 0000000000..d874b4906a
Binary files /dev/null and b/docs/img/plotmeshtally.png differ
diff --git a/docs/source/devguide/index.rst b/docs/source/devguide/index.rst
index 98df54b17e..1eb7436d03 100644
--- a/docs/source/devguide/index.rst
+++ b/docs/source/devguide/index.rst
@@ -17,3 +17,4 @@ as debugging.
workflow
xml-fortran
statepoint
+ voxel
diff --git a/docs/source/devguide/voxel.rst b/docs/source/devguide/voxel.rst
new file mode 100644
index 0000000000..8ce0c6ee87
--- /dev/null
+++ b/docs/source/devguide/voxel.rst
@@ -0,0 +1,54 @@
+.. _devguide_voxel:
+
+=====================================
+Voxel Plot Binary File Specifications
+=====================================
+
+----------
+Revision 1
+----------
+
+**integer(4) n_voxels_x**
+
+ Number of voxels in the x direction
+
+**integer(4) n_voxels_y**
+
+ Number of voxels in the y direction
+
+**integer(4) n_voxels_z**
+
+ Number of voxels in the z direction
+
+**real(8) width_voxel_x**
+
+ Width of voxels in the x direction
+
+**real(8) width_voxel_y**
+
+ Width of voxels in the y direction
+
+**real(8) width_voxel_z**
+
+ Width of voxels in the z direction
+
+**real(8) lower_left_x**
+
+ Lower left x point of the voxel grid
+
+**real(8) lower_left_y**
+
+ Lower left y point of the voxel grid
+
+**real(8) lower_left_z**
+
+ Lower left z point of the voxel grid
+
+*do x = 1, n_voxels_x*
+ *do y = 1, n_voxels_y*
+ *do z = 1, n_voxels_z*
+
+ **integer(4) id**
+
+ Cell or material id number at this voxel center. Set to -1 when
+ cell not_found.
diff --git a/docs/source/usersguide/beginners.rst b/docs/source/usersguide/beginners.rst
index 9a791ac521..69f44c2ca1 100644
--- a/docs/source/usersguide/beginners.rst
+++ b/docs/source/usersguide/beginners.rst
@@ -109,7 +109,9 @@ should be comfortable working in a command line environment. There are many
resources online for learning command line environments. If you are using Linux
or Mac OS X (also Unix-derived), `this tutorial
`_ will help you get acquainted with
-commonly-used commands.
+commonly-used commands. It is also helpful to be familiar with `Python
+`_, as most of the post-processing utilities provided
+with OpenMC rely on it for data manipulation and results visualization.
OpenMC uses a version control software called `git`_ to keep track of changes to
the code, document bugs and issues, and other development tasks. While you don't
diff --git a/docs/source/usersguide/index.rst b/docs/source/usersguide/index.rst
index a3b1583653..675ed4081a 100644
--- a/docs/source/usersguide/index.rst
+++ b/docs/source/usersguide/index.rst
@@ -14,4 +14,5 @@ essential aspects of using OpenMC to perform neutronic simulations.
beginners
install
input
+ processing
troubleshoot
diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst
index af8d94e71d..3571f500f1 100644
--- a/docs/source/usersguide/input.rst
+++ b/docs/source/usersguide/input.rst
@@ -314,7 +314,7 @@ attributes/sub-elements:
*Default*: watt
:parameters:
- For a "monoenergetic" energy distribution, ``parameters`` should not be
+ For a "monoenergetic" energy distribution, ``parameters`` should be
given as the energy in MeV of the source sites.
For a "watt" energy distribution, ``parameters`` should be given as two
@@ -941,19 +941,27 @@ tallies. This element should be followed by "true" or "false".
*Default*: false
+.. _usersguide_plotting:
+
--------------------------------------------
Geometry Plotting Specification -- plots.xml
--------------------------------------------
-A basic 2D plotting capability is available in OpenMC by creating a plots.xml
+Basic plotting capabilities are available in OpenMC by creating a plots.xml
file and subsequently running with the command-line flag ``-plot``. The root
-element of the plots.xml is simply ```` and any number output figures can
-be defined with ```` sub-elements.
+element of the plots.xml is simply ```` and any number output plots can
+be defined with ```` sub-elements. Two plot types are currently
+implemented in openMC:
+
+* ``slice`` 2D pixel plot along one of the major axes. Produces a PPM image file.
+* ``voxel`` 3D voxel data dump. Produces a binary file containing voxel xyz position and cell or material id.
+
```` Element
------------------
-Each plot must contain a combination of the following attributes or sub-elements:
+Each plot must contain a combination of the following attributes or
+sub-elements:
:id:
The unique ``id`` of the plot.
@@ -967,7 +975,9 @@ Each plot must contain a combination of the following attributes or sub-elements
:color:
Keyword for plot coloring. This can only be either ``cell`` or ``mat``,
- which colors regions by cells and materials, respectively.
+ which colors regions by cells and materials, respectively. For voxel plots,
+ this determines which id (cell or material) is associated with each
+ position.
*Default*: ``cell``
@@ -985,60 +995,75 @@ Each plot must contain a combination of the following attributes or sub-elements
*Default*: None - Required entry
:type:
- Keyword for type of plot to be produced. Currently only ``slice`` plots are
- implemented, which create 2D pixel maps saved in the PPM file format. PPM
- files can be displayed in most viewers (e.g. the default Gnome viewer,
- IrfanView, etc.).
+ Keyword for type of plot to be produced. Currently only "slice" and
+ "voxel" plots are implemented. The "slice" plot type creates 2D pixel
+ maps saved in the PPM file format. PPM files can be displayed in most
+ viewers (e.g. the default Gnome viewer, IrfanView, etc.). The "voxel"
+ plot type produces a binary datafile containing voxel grid positioning and
+ the cell or material (specified by the ``color`` tag) at the center of each
+ voxel. These datafiles can be processed into 3D SILO files using the
+ ``voxel.py`` utility provided with the OpenMC source, and subsequently
+ viewed with a 3D viewer such as VISIT or Paraview. See the
+ :ref:`devguide_voxel` for information about the datafile structure.
.. note:: Since the PPM format is saved without any kind of compression,
the resulting file sizes can be quite large. Saving the image in
the PNG format can often times reduce the file size by orders of
- magnitude without any loss of image quality.
+ magnitude without any loss of image quality. Likewise,
+ high-resolution voxel files produced by OpenMC can be quite large,
+ but the equivalent SILO files will by significantly smaller.
*Default*: "slice"
-```` elements of ``type`` "slice" also contain the following attributes or
-sub-elements:
+```` elements of ``type`` "slice" and "voxel" must contain the ``pixels``
+attribute or sub-element:
+
+ :pixels:
+ Specifies the number of pixes or voxels to be used along each of the basis
+ directions for "slice" and "voxel" plots, respectively. Should be two or
+ three integers separated by spaces.
+
+ .. warning:: The ``pixels`` input determines the output file size. For the
+ PPM format, 10 million pixels will result in a file just under
+ 30 MB in size. A 10 million voxel binary file will be around
+ 40 MB.
+
+ .. warning:: If the aspect ratio defined in ``pixels`` does not match the
+ aspect ratio defined in ``width`` the plot may appear stretched
+ or squeezed.
+
+ .. warning:: Geometry features along a basis direction smaller than
+ ``width``/``pixels`` along that basis direction may not appear
+ in the plot.
+
+ *Default*: None - Required entry for "slice" and "voxel" plots
+
+```` elements of ``type`` "slice" can also contain the following
+attributes or sub-elements. These are not used in "voxel" plots:
:basis:
- Keyword specifying the plane of the plot for ``slice`` type plots. Can be
+ Keyword specifying the plane of the plot for "slice" type plots. Can be
one of: "xy", "xz", "yz".
*Default*: "xy"
- :pixels:
- Specifies the number of pixes to be used along each of the basis directions
- for "slice" plots. Should be two integers separated by spaces.
-
- .. warning:: The ``pixels`` input determines the output file size. For the PPM
- format, 10 million pixels will result in a file just under 30 MB in
- size.
-
- .. warning:: If the aspect ratio defined in ``pixels`` does not match the aspect
- ratio defined in ``width`` the plot may appear stretched or squeezed.
-
- .. warning:: Geometry features along a basis direction smaller than ``width``/``pixels``
- along that basis direction may not appear in the plot.
-
- *Default*: None - Required entry for "slice" plots
-
:background:
- Specifies the RGB color of the regions where no OpenMC cell can be found. Should
- be three integers separated by spaces.
+ Specifies the RGB color of the regions where no OpenMC cell can be found.
+ Should be three integers separated by spaces.
*Default*: 0 0 0 (white)
:col_spec:
- Any number of this optional tag may be included in each ```` element, which can
- override the default random colors for cells or materials. Each ``col_spec``
- element must contain ``id`` and ``rgb`` sub-elements.
+ Any number of this optional tag may be included in each ```` element,
+ which can override the default random colors for cells or materials. Each
+ ``col_spec`` element must contain ``id`` and ``rgb`` sub-elements.
:id:
Specifies the cell or material unique id for the color specification.
:rgb:
- Specifies the custom color for the cell or material. Should be 3 integers separated
- by spaces.
+ Specifies the custom color for the cell or material. Should be 3 integers
+ separated by spaces.
As an example, if your plot is colored by material and you want material 23
to be blue, the corresponding ``col_spec`` element would look like:
@@ -1051,17 +1076,18 @@ sub-elements:
:mask:
The special ``mask`` sub-element allows for the selective plotting of *only*
- user-specified cells or materials. Only one ``mask`` element is allowed per ``plot``
- element, and it must contain as attributes or sub-elements a background masking color and
- a list of cells or materials to plot:
+ user-specified cells or materials. Only one ``mask`` element is allowed per
+ ``plot`` element, and it must contain as attributes or sub-elements a
+ background masking color and a list of cells or materials to plot:
:components:
- List of unique ``id`` numbers of the cells or materials to plot. Should be any number
- of integers separated by spaces.
+ List of unique ``id`` numbers of the cells or materials to plot. Should be
+ any number of integers separated by spaces.
:background:
- Color to apply to all cells or materials not in the ``components`` list of cells or
- materials to plot. This overrides any ``col_spec`` color specifications.
+ Color to apply to all cells or materials not in the ``components`` list of
+ cells or materials to plot. This overrides any ``col_spec`` color
+ specifications.
*Default*: None
diff --git a/docs/source/usersguide/processing.rst b/docs/source/usersguide/processing.rst
new file mode 100644
index 0000000000..4e9daa4072
--- /dev/null
+++ b/docs/source/usersguide/processing.rst
@@ -0,0 +1,360 @@
+.. _usersguide_processing:
+
+=================================
+Data Processing and Visualization
+=================================
+
+This section is intended to explain in detail the recommended procedures for
+carrying out common tasks with OpenMC. While several utilities of varying
+complexity are provided to help automate the process, in many cases it will be
+extremely beneficial to do some coding in Python to quickly obtain results. In
+these cases, and for many of the provided utilities, it is necessary for your
+Python installation to contain:
+
+* [1]_ `Numpy `_
+* [1]_ `Scipy `_
+* [2]_ `h5py `_
+* [3]_ `Matplotlib `_
+* [3]_ `Silomesh `_
+* [3]_ `VTK `_
+* [4]_ `PyQt `_
+
+Most of these are easily obtainable in Ubuntu through the package manager, or
+are easily installed with distutils.
+
+.. [1] Required for tally data extraction from statepoints with statepoint.py
+.. [2] Required only if reading HDF5 statepoint files.
+.. [3] Optional for plotting utilities
+.. [4] Optional for interactive GUIs
+
+----------------------
+Geometry Visualization
+----------------------
+
+Geometry plotting is carried out by creating a plots.xml, specifying plots, and
+running OpenMC with the -plot or -p command-line option (See
+:ref:`usersguide_plotting`).
+
+Plotting in 2D
+--------------
+
+.. image:: ../../img/atr.png
+ :height: 200px
+
+After running OpenMC to obtain PPM files, images should be saved to another
+format before using them elsewhere. This cuts down the size of the file by
+orders of magnitude. Most image viewers and editors that can view PPM images
+can also save to other formats (e.g. `Gimp `_, `IrfanView
+`_, etc.). However, more likey the user will want to
+convert to another format on the command line. This is easily accomplished with
+the ``convert`` command available on most linux distributions as part of the
+`ImageMagick `_ package. (On
+Ubuntu: ``sudo apt-get install imagemagick``). Images are then converted like:
+
+.. code-block:: sh
+
+ convert plot.ppm plot.png
+
+Plotting in 3D
+--------------
+
+.. image:: ../../img/3dgeomplot.png
+ :height: 200px
+
+The binary VOXEL files output by OpenMC can not be viewed directly by any
+existing viewers. In order to view them, they must be converted into a standard
+mesh format that can be viewed in ParaView, Visit, etc. The provided utility
+voxel.py accomplishes this for SILO:
+
+.. code-block:: sh
+
+ /src/utils/voxel.py myplot.voxel -o output.silo
+
+and VTK file formats:
+
+.. code-block:: sh
+
+ /src/utils/myplot.voxel --vtk -o output.vti
+
+To use this utility you need either
+
+* `Silomesh `_
+
+or
+
+* `VTK `_ with python bindings - On Ubuntu, these are easily obtained with ``sudo apt-get install python-vtk``
+
+Users can process the binary into any other format if desired by following the
+example of voxel.py. For the binary file structure, see :ref:`devguide_voxel`.
+
+.. note:: 3D voxel plotting can be very computer intensive for the viewing
+ program (Visit, Paraview, etc.) if the number of voxels is large
+ (>10million or so). Thus if you want an accurate picture that
+ renders smoothly, consider using only one voxel in a certain
+ direction. For instance, the 3D pin lattice figure above was generated
+ with a 500x500x1 voxel mesh, which allows for resolution of the
+ cylinders without wasting too many voxels on the axial dimension.
+
+
+-------------------
+Tally Visualization
+-------------------
+
+Tally results are saved in both a text file (tallies.out) as well as a binary
+statepoint file. While the tallies.out file may be fine for simple tallies, in
+many cases the user requires more information about the tally or the run, or
+has to deal with a large number of result values (e.g. for mesh tallies). In
+these cases, extracting data from the statepoint file via Python scripting is
+the preferred method of data analysis and visualization.
+
+Data Extraction
+---------------
+
+A great deal of information is available in statepoint files (See
+:ref:`devguide_statepoint`), most of which is easily extracted by the provided
+utility statepoint.py. This utility provides a Python class to load statepoints
+and extract data - it is used in many of the provided plotting utilities, and
+can be used in user-created scripts to carry out manipulations of the data. To
+read tallies using this utility, make sure statepoint.py is in your PYTHONPATH,
+and then import the class, instantiate it, and call read_results:
+
+.. code-block:: python
+
+ from statepoint import StatePoint
+ sp = StatePoint('statepoint.100.binary')
+ sp.read_results()
+
+At this point the user can extract entire scores from tallies into a data
+dictionary containing numpy arrays:
+
+.. code-block:: python
+
+ tallyid = 1
+ score = 'flux'
+ data = sp.extract_results(tallyid, score)
+ means = data['means']
+ print data.keys()
+
+The results from this function contain all filter bins (all mesh points, all
+energy groups, etc.), which can be reshaped with the bin ordering also contained
+in the output dictionary. This is the best choice of output for easily
+integrating ranges of data.
+
+Alternatively the user can extract specific values for a single score/filter
+combination:
+
+.. code-block:: python
+
+ tallyid = 1
+ score = 'flux'
+ filters = [('mesh', (1, 1, 5)), ('energyin', 0)]
+ value, error = sp.get_value(tallyid, filters, score)
+
+In the future more documentaion may become available here for statepoint.py and
+the data extraction functions of StatePoint objects. However, for now it is up
+to the user to explore the classes in statepoint.py to discover what data is
+available in StatePoint objects (we highly recommend interactively exploring
+with `IPython `_). Many exmaples can be found by looking
+through the other utilies that use statepoint.py, and a few common visualization
+tasks will be described here in the following sections.
+
+Plotting in 2D
+--------------
+
+.. image:: ../../img/plotmeshtally.png
+ :height: 200px
+
+For simple viewing of 2D slices of a mesh plot, the utility plot_mesh_tally.py
+is provided. This utility provides an interactive GUI to explore and plot
+mesh tallies for any scores and filter bins. It requires statepoint.py, as well
+as `PyQt `_.
+
+.. image:: ../../img/fluxplot.png
+ :height: 200px
+
+Alternatively, the user can write their own Python script to manipulate the data
+appropriately. Consider a run where the first tally contains a 105x105x20 mesh
+over a small core, with a flux score and two energyin filter bins. To explicitly
+extract the data and create a plot with gnuplot, the following script can be
+used. The script operates in several steps for clarity, and is not necessarily
+the most efficient way to extract data from large mesh tallies. This creates the
+two heatmaps in the previous figure.
+
+.. code-block:: python
+
+ #!/usr/bin/env python
+
+ import os
+
+ import statepoint
+
+ # load and parse the statepoint file
+ sp = statepoint.StatePoint('statepoint.300.binary')
+ sp.read_results()
+
+ tallyid = 0 # This is tally 1
+ score = 0 # This corresponds to flux (see tally.scores)
+
+ # get mesh dimensions
+ meshid = sp.tallies[tallyid].filters['mesh'].bins[0]
+ for i,m in enumerate(sp.meshes):
+ if m.id == meshid:
+ mesh = m
+ break
+ nx,ny,nz = mesh.dimension
+
+ # loop through mesh and extract values to python dictionaries
+ thermal = {}
+ fast = {}
+ for x in range(1,nx+1):
+ for y in range(1,ny+1):
+ for z in range(1,nz+1):
+ val,err = sp.get_value(tallyid,
+ [('mesh',(x,y,z)),('energyin',0)],
+ score)
+ thermal[(x,y,z)] = val
+ val,err = sp.get_value(tallyid,
+ [('mesh',(x,y,z)),('energyin',1)],
+ score)
+ fast[(x,y,z)] = val
+
+ # sum up the axial values and write datafile for gnuplot
+ with open('meshdata.dat','w') as fh:
+ for x in range(1,nx+1):
+ for y in range(1,ny+1):
+ thermalval = 0.
+ fastval = 0.
+ for z in range(1,nz+1):
+ thermalval += thermal[(x,y,z)]
+ fastval += fast[(x,y,z)]
+ fh.write("{} {} {} {}\n".format(x,y,thermalval,fastval))
+
+ # write gnuplot file
+ with open('tmp.gnuplot','w') as fh:
+ fh.write(r"""set terminal png size 1000 400
+ set output 'fluxplot.png'
+ set nokey
+ set autoscale fix
+ set multiplot layout 1,2 title "Pin Mesh Flux Tally"
+ set title "Thermal"
+ plot 'meshdata.dat' using 1:2:3 with image
+ set title "Fast"
+ plot 'meshdata.dat' using 1:2:4 with image
+ """)
+
+ # make plot
+ os.system("gnuplot < tmp.gnuplot")
+
+Plotting in 3D
+--------------
+
+.. image:: ../../img/3dcore.png
+ :height: 200px
+
+As with 3D plots of the geometry, meshtally data needs to be put into a standard
+format for viewing. The utility statepoint_3d.py is provided to accomplish this
+for both VTK and SILO. By default statepoint_3d.py processes a statepoint into a
+3D file with all mesh tallies and filter/score combinations,
+
+.. code-block:: sh
+
+ /src/utils/statepoint_3d.py -o output.silo
+ /src/utils/statepoint_3d.py --vtk -o output.vtm
+
+but it also provides several command-line options to selectively process only
+certain data arrays in order to keep file sizes down.
+
+.. code-block:: sh
+
+ /src/utils/statepoint_3d.py --tallies 2,4 --scores 4.1,4.3 -o output.silo
+ /src/utils/statepoint_3d.py --filters 2.energyin.1 --vtk -o output.vtm
+
+All available options for specifying a subset of tallies, scores, and filters
+can be listed with the ``--list`` or ``-l`` command line options.
+
+.. note:: Note that while SILO files can contain multiple meshes in one file,
+ VTK needs to use a multi-block dataset, which stores each mesh piece
+ in a different file in a subfolder. All meshes can be loaded at once
+ with the main VTM file, or each VTI file in the subfolder can be
+ loaded individually.
+
+Alternatively, the user can write their own Python script to manipulate the data
+appropriately before insertion into a SILO or VTK file. For instance, if the
+data has been extracted as was done in the 2D plotting example script above, a
+SILO file can be created with:
+
+.. code-block:: python
+
+ import silomesh as sm
+ sm.init_silo("fluxtally.silo")
+ sm.init_mesh('tally_mesh',*mesh.dimension, *mesh.lower_left, *mesh.upper_right)
+ sm.init_var('flux_tally_thermal')
+ for x in range(1,nx+1):
+ for y in range(1,ny+1):
+ for z in range(1,nz+1):
+ sm.set_value(float(thermal[(x,y,z)]),x,y,z)
+ sm.finalize_var()
+ sm.init_var('flux_tally_fast')
+ for x in range(1,nx+1):
+ for y in range(1,ny+1):
+ for z in range(1,nz+1):
+ sm.set_value(float(fast[(x,y,z)]),x,y,z)
+ sm.finalize_var()
+ sm.finalize_mesh()
+ sm.finalize_silo()
+
+and the equivalent VTK file with:
+
+.. code-block:: python
+
+ import vtk
+
+ grid = vtk.vtkImageData()
+ grid.SetDimensions(nx+1,ny+1,nz+1)
+ grid.SetOrigin(*mesh.lower_left)
+ grid.SetSpacing(*mesh.width)
+
+ # vtk cell arrays have x on the inners, so we need to reorder the data
+ idata = {}
+ for x in range(nx):
+ for y in range(ny):
+ for z in range(nz):
+ i = z*nx*ny + y*nx + x
+ idata[i] = (x,y,z)
+
+ vtkfastdata = vtk.vtkDoubleArray()
+ vtkfastdata.SetName("fast")
+ for i in range(nx*ny*nz):
+ vtkfastdata.InsertNextValue(fast[idata[i]])
+
+ vtkthermaldata = vtk.vtkDoubleArray()
+ vtkthermaldata.SetName("thermal")
+ for i in range(nx*ny*nz):
+ vtkthermaldata.InsertNextValue(thermal[idata[i]])
+
+ grid.GetCellData().AddArray(vtkfastdata)
+ grid.GetCellData().AddArray(vtkthermaldata)
+
+ writer = vtk.vtkXMLImageDataWriter()
+ writer.SetInput(grid)
+ writer.SetFileName('tally.vti')
+ writer.Write()
+
+Getting Data into MATLAB
+------------------------
+
+There is currently no front-end utility to dump tally data to MATLAB files, but
+the process is straightforward. First extract the data using a custom Python
+script with statepoint.py, put the data into appropriately-shaped numpy arrays,
+and then use the `Scipy MATLAB IO routines
+`_ to save to a MAT
+file. Note that the data contained in the output from
+``StatePoint.extract_result`` is already in a Numpy array that can be reshaped
+and dumped to MATLAB in one step.
+
+
+
+
+
+
+
diff --git a/src/constants.F90 b/src/constants.F90
index 9ff07788e8..04f757db9c 100644
--- a/src/constants.F90
+++ b/src/constants.F90
@@ -48,13 +48,16 @@ module constants
! ============================================================================
! PHYSICAL CONSTANTS
+ ! Values here are from the Committee on Data for Science and Technology
+ ! (CODATA) 2010 recommendation (doi:10.1103/RevModPhys.84.1527).
+
real(8), parameter :: &
PI = 3.1415926535898_8, & ! pi
- MASS_NEUTRON = 1.0086649156, & ! mass of a neutron
- MASS_PROTON = 1.00727646677, & ! mass of a proton
- AMU = 1.66053873e-27, & ! 1 amu in kg
- N_AVOGADRO = 0.602214179, & ! Avogadro's number in 10^24/mol
- K_BOLTZMANN = 8.617342e-11, & ! Boltzmann constant in MeV/K
+ MASS_NEUTRON = 1.008664916, & ! mass of a neutron in amu
+ MASS_PROTON = 1.007276466812, & ! mass of a proton in amu
+ AMU = 1.660538921e-27, & ! 1 amu in kg
+ N_AVOGADRO = 0.602214129, & ! Avogadro's number in 10^24/mol
+ K_BOLTZMANN = 8.6173324e-11, & ! Boltzmann constant in MeV/K
INFINITY = huge(0.0_8), & ! positive infinity
ZERO = 0.0_8, &
ONE = 1.0_8, &
diff --git a/src/geometry.F90 b/src/geometry.F90
index 6859712886..ecf1698ff9 100644
--- a/src/geometry.F90
+++ b/src/geometry.F90
@@ -1114,7 +1114,7 @@ contains
! Check is calculated distance is new minimum
if (d < dist) then
- if (abs(d - dist)/dist >= FP_PRECISION) then
+ if (abs(d - dist)/dist >= FP_REL_PRECISION) then
dist = d
surface_crossed = -cl % surfaces(i)
lattice_crossed = NONE
@@ -1156,7 +1156,8 @@ contains
! point precision.
if (d < dist) then
- if (abs(d - dist)/dist >= FP_REL_PRECISION) then
+ if (abs(d - dist)/dist >= FP_REL_PRECISION &
+ .and. abs(d - dist) >= FP_PRECISION) then
dist = d
if (u > 0) then
lattice_crossed = LATTICE_RIGHT
@@ -1177,7 +1178,8 @@ contains
end if
if (d < dist) then
- if (abs(d - dist)/dist >= FP_REL_PRECISION) then
+ if (abs(d - dist)/dist >= FP_REL_PRECISION &
+ .and. abs(d - dist) >= FP_PRECISION) then
dist = d
if (v > 0) then
lattice_crossed = LATTICE_FRONT
@@ -1201,7 +1203,8 @@ contains
end if
if (d < dist) then
- if (abs(d - dist)/dist >= FP_REL_PRECISION) then
+ if (abs(d - dist)/dist >= FP_REL_PRECISION &
+ .and. abs(d - dist) >= FP_PRECISION) then
dist = d
if (w > 0) then
lattice_crossed = LATTICE_TOP
diff --git a/src/global.F90 b/src/global.F90
index 5042ac87ec..0327e0fd19 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -10,7 +10,7 @@ module global
use material_header, only: Material
use mesh_header, only: StructuredMesh
use particle_header, only: Particle
- use plot_header, only: PlotSlice
+ use plot_header, only: ObjectPlot
use set_header, only: SetInt
use source_header, only: ExtSource
use tally_header, only: TallyObject, TallyMap, TallyResult
@@ -41,7 +41,7 @@ module global
type(Lattice), allocatable, target :: lattices(:)
type(Surface), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
- type(PlotSlice), allocatable, target :: plots(:)
+ type(ObjectPlot),allocatable, target :: plots(:)
! Size of main arrays
integer :: n_cells ! # of cells
diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90
index 29c36aaae0..8189a56edf 100644
--- a/src/hdf5_interface.F90
+++ b/src/hdf5_interface.F90
@@ -1195,6 +1195,7 @@ contains
dims(1) = restart_batch
call h5ltread_dataset_double_f(hdf5_state_point, "k_batch", &
k_batch(1:restart_batch), dims, hdf5_err)
+ dims(1) = restart_batch*gen_per_batch
call h5ltread_dataset_double_f(hdf5_state_point, "entropy", &
entropy(1:restart_batch*gen_per_batch), dims, hdf5_err)
call hdf5_read_double(hdf5_state_point, "k_col_abs", k_col_abs)
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 8a927ef9c9..c07fe6a5c3 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -2135,7 +2135,7 @@ contains
integer n_cols, col_id
logical :: file_exists ! does plots.xml file exist?
character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml
- type(PlotSlice), pointer :: pl => null()
+ type(ObjectPlot), pointer :: pl => null()
! Check if plots.xml exists
filename = trim(path_input) // "plots.xml"
@@ -2169,20 +2169,55 @@ contains
call fatal_error()
end if
- ! Set output file path
- pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // &
- "_" // trim(plot_(i) % filename) // ".ppm"
-
- ! Copy plot pixel size
- if (size(plot_(i) % pixels) == 2) then
- pl % pixels = plot_(i) % pixels
- else
- message = " must be length 2 in plot " // to_str(pl % id)
+ ! Copy plot type
+ select case (plot_(i) % type)
+ case ("slice")
+ pl % type = PLOT_TYPE_SLICE
+ case ("voxel")
+ pl % type = PLOT_TYPE_VOXEL
+ case default
+ message = "Unsupported plot type '" // trim(plot_(i) % type) &
+ // "' in plot " // trim(to_str(pl % id))
call fatal_error()
+ end select
+
+ ! Set output file path
+ select case (pl % type)
+ case (PLOT_TYPE_SLICE)
+ pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // &
+ "_" // trim(plot_(i) % filename) // ".ppm"
+ case (PLOT_TYPE_VOXEL)
+ pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // &
+ "_" // trim(plot_(i) % filename) // ".voxel"
+ end select
+
+ ! Copy plot pixel size
+ if (pl % type == PLOT_TYPE_SLICE) then
+ if (size(plot_(i) % pixels) == 2) then
+ pl % pixels(1) = plot_(i) % pixels(1)
+ pl % pixels(2) = plot_(i) % pixels(2)
+ else
+ message = " must be length 2 in slice plot " // &
+ trim(to_str(pl % id))
+ call fatal_error()
+ end if
+ else if (pl % type == PLOT_TYPE_VOXEL) then
+ if (size(plot_(i) % pixels) == 3) then
+ pl % pixels = plot_(i) % pixels
+ else
+ message = " must be length 3 in voxel plot " // &
+ trim(to_str(pl % id))
+ call fatal_error()
+ end if
end if
! Copy plot background color
if (associated(plot_(i) % background)) then
+ if (pl % type == PLOT_TYPE_VOXEL) then
+ message = "Background color ignored in voxel plot " // &
+ trim(to_str(pl % id))
+ call warning()
+ end if
if (size(plot_(i) % background) == 3) then
pl % not_found % rgb = plot_(i) % background
else
@@ -2194,30 +2229,22 @@ contains
pl % not_found % rgb = (/ 255, 255, 255 /)
end if
- ! Copy plot type
- select case (plot_(i) % type)
- case ("slice")
- pl % type = PLOT_TYPE_SLICE
- case default
- message = "Unsupported plot type '" // plot_(i) % type &
- // "' in plot " // trim(to_str(pl % id))
- call fatal_error()
- end select
-
! Copy plot basis
- select case (plot_(i) % basis)
- case ("xy")
- pl % basis = PLOT_BASIS_XY
- case ("xz")
- pl % basis = PLOT_BASIS_XZ
- case ("yz")
- pl % basis = PLOT_BASIS_YZ
- case default
- message = "Unsupported plot basis '" // plot_(i) % basis &
- // "' in plot " // trim(to_str(pl % id))
- call fatal_error()
- end select
-
+ if (pl % type == PLOT_TYPE_SLICE) then
+ select case (plot_(i) % basis)
+ case ("xy")
+ pl % basis = PLOT_BASIS_XY
+ case ("xz")
+ pl % basis = PLOT_BASIS_XZ
+ case ("yz")
+ pl % basis = PLOT_BASIS_YZ
+ case default
+ message = "Unsupported plot basis '" // plot_(i) % basis &
+ // "' in plot " // trim(to_str(pl % id))
+ call fatal_error()
+ end select
+ end if
+
! Copy plotting origin
if (size(plot_(i) % origin) == 3) then
pl % origin = plot_(i) % origin
@@ -2228,15 +2255,23 @@ contains
end if
! Copy plotting width
- if (size(plot_(i) % width) == 3) then
- pl % width = plot_(i) % width
- else if (size(plot_(i) % width) == 2) then
- pl % width(1) = plot_(i) % width(1)
- pl % width(2) = plot_(i) % width(2)
- else
- message = "Bad plot width " &
- // "in plot " // trim(to_str(pl % id))
- call fatal_error()
+ if (pl % type == PLOT_TYPE_SLICE) then
+ if (size(plot_(i) % width) == 2) then
+ pl % width(1) = plot_(i) % width(1)
+ pl % width(2) = plot_(i) % width(2)
+ else
+ message = " must be length 2 in slice plot " // &
+ trim(to_str(pl % id))
+ call fatal_error()
+ end if
+ else if (pl % type == PLOT_TYPE_VOXEL) then
+ if (size(plot_(i) % width) == 3) then
+ pl % width = plot_(i) % width
+ else
+ message = " must be length 3 in voxel plot " // &
+ trim(to_str(pl % id))
+ call fatal_error()
+ end if
end if
! Copy plot color type and initialize all colors randomly
@@ -2269,6 +2304,13 @@ contains
! Copy user specified colors
if (associated(plot_(i) % col_spec_)) then
+
+ if (pl % type == PLOT_TYPE_VOXEL) then
+ message = "Color specifications ignored in voxel plot " // &
+ trim(to_str(pl % id))
+ call warning()
+ end if
+
n_cols = size(plot_(i) % col_spec_)
do j = 1, n_cols
if (size(plot_(i) % col_spec_(j) % rgb) /= 3) then
@@ -2308,6 +2350,12 @@ contains
! Deal with masks
if (associated(plot_(i) % mask_)) then
+ if (pl % type == PLOT_TYPE_VOXEL) then
+ message = "Mask ignored in voxel plot " // &
+ trim(to_str(pl % id))
+ call warning()
+ end if
+
select case(size(plot_(i) % mask_))
case default
message = "Mutliple masks" // &
diff --git a/src/output.F90 b/src/output.F90
index 9915a27474..edb83be0da 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -1308,7 +1308,7 @@ contains
subroutine print_plot()
integer :: i ! loop index for plots
- type(PlotSlice), pointer :: pl => null()
+ type(ObjectPlot), pointer :: pl => null()
! Display header for plotting
call header("PLOTTING SUMMARY")
@@ -1316,20 +1316,34 @@ contains
do i = 1, n_plots
pl => plots(i)
- ! Write plot id
+ ! Plot id
write(ou,100) "Plot ID:", trim(to_str(pl % id))
+
+ ! Plot type
+ if (pl % type == PLOT_TYPE_SLICE) then
+ write(ou,100) "Plot Type:", "Slice"
+ else if (pl % type == PLOT_TYPE_VOXEL) then
+ write(ou,100) "Plot Type:", "Voxel"
+ end if
- ! Write plotting origin
+ ! Plot parameters
write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // &
" " // trim(to_str(pl % origin(2))) // " " // &
trim(to_str(pl % origin(3)))
-
- ! Write plotting width
if (pl % type == PLOT_TYPE_SLICE) then
-
write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
" " // trim(to_str(pl % width(2)))
- write(ou,100) "Coloring:", trim(to_str(pl % color_by))
+ else if (pl % type == PLOT_TYPE_VOXEL) then
+ write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
+ " " // trim(to_str(pl % width(2))) // &
+ " " // trim(to_str(pl % width(3)))
+ end if
+ if (pl % color_by == PLOT_COLOR_CELLS) then
+ write(ou,100) "Coloring:", "Cells"
+ else if (pl % color_by == PLOT_COLOR_MATS) then
+ write(ou,100) "Coloring:", "Materials"
+ end if
+ if (pl % type == PLOT_TYPE_SLICE) then
select case (pl % basis)
case (PLOT_BASIS_XY)
write(ou,100) "Basis:", "xy"
@@ -1340,6 +1354,9 @@ contains
end select
write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // &
trim(to_str(pl % pixels(2)))
+ else if (pl % type == PLOT_TYPE_VOXEL) then
+ write(ou,100) "Voxels:", trim(to_str(pl % pixels(1))) // " " // &
+ trim(to_str(pl % pixels(2))) // " " // trim(to_str(pl % pixels(3)))
end if
write(ou,*)
diff --git a/src/physics.F90 b/src/physics.F90
index 76c4141cdd..3000a2dbce 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -790,7 +790,7 @@ contains
kT = nuc % kT
! Check if energy is above threshold
- if (p % E >= FREE_GAS_THRESHOLD * kT) then
+ if (p % E >= FREE_GAS_THRESHOLD * kT .and. nuc % awr > ONE) then
v_target = ZERO
return
end if
diff --git a/src/plot.F90 b/src/plot.F90
index 85fc2c8ac6..0b68f4c219 100644
--- a/src/plot.F90
+++ b/src/plot.F90
@@ -24,7 +24,7 @@ contains
subroutine run_plot()
integer :: i ! loop index for plots
- type(PlotSlice), pointer :: pl => null()
+ type(ObjectPlot), pointer :: pl => null()
do i = 1, n_plots
pl => plots(i)
@@ -36,11 +36,60 @@ contains
if (pl % type == PLOT_TYPE_SLICE) then
! create 2d image
call create_ppm(pl)
+ else if (pl % type == PLOT_TYPE_VOXEL) then
+ ! create dump for 3D silomesh utility script
+ call create_3d_dump(pl)
end if
end do
end subroutine run_plot
+!===============================================================================
+! POSITION_RGB computes the red/green/blue values for a given plot with the
+! current particle's position
+!===============================================================================
+
+ subroutine position_rgb(pl, rgb, id)
+
+ type(ObjectPlot), pointer, intent(in) :: pl
+ integer, intent(out) :: rgb(3)
+ integer, intent(out) :: id
+
+ logical :: found_cell
+ type(Cell), pointer :: c => null()
+
+ call deallocate_coord(p % coord0 % next)
+ p % coord => p % coord0
+
+ call find_cell(found_cell)
+
+ if (.not. found_cell) then
+ ! If no cell, revert to default color
+ rgb = pl % not_found % rgb
+ id = -1
+ else
+ if (pl % color_by == PLOT_COLOR_MATS) then
+ ! Assign color based on material
+ c => cells(p % coord % cell)
+ id = materials(c % material) % id
+ if (c % material == MATERIAL_VOID) then
+ ! By default, color void cells white
+ rgb = 255
+ else
+ rgb = pl % colors(c % material) % rgb
+ end if
+ else if (pl % color_by == PLOT_COLOR_CELLS) then
+ ! Assign color based on cell
+ rgb = pl % colors(p % coord % cell) % rgb
+ id = cells(p % coord % cell) % id
+ else
+ rgb = 0
+ id = -1
+ end if
+ end if
+
+ end subroutine position_rgb
+
!===============================================================================
! CREATE_PPM creates an image based on user input from a plots.xml
! specification in the portable pixmap format (PPM)
@@ -48,18 +97,17 @@ contains
subroutine create_ppm(pl)
- type(PlotSlice), pointer :: pl
+ type(ObjectPlot), pointer :: pl
integer :: in_i
integer :: out_i
integer :: x, y ! pixel location
- integer :: r, g, b ! colors (red, green, blue) from 0-255
+ integer :: rgb(3) ! colors (red, green, blue) from 0-255
+ integer :: id
real(8) :: in_pixel
real(8) :: out_pixel
real(8) :: xyz(3)
- logical :: found_cell
type(Image) :: img
- type(Cell), pointer :: c => null()
! Initialize and allocate space for image
call init_image(img)
@@ -101,44 +149,11 @@ contains
do y = 1, img % height
do x = 1, img % width
- call deallocate_coord(p % coord0 % next)
- p % coord => p % coord0
-
- call find_cell(found_cell)
-
- if (.not. found_cell) then
- ! If no cell, revert to default color
- r = pl % not_found % rgb(1)
- g = pl % not_found % rgb(2)
- b = pl % not_found % rgb(3)
- else
- if (pl % color_by == PLOT_COLOR_MATS) then
- ! Assign color based on material
- c => cells(p % coord % cell)
- if (c % material == MATERIAL_VOID) then
- ! By default, color void cells white
- r = 255
- g = 255
- b = 255
- else
- r = pl % colors(c % material) % rgb(1)
- g = pl % colors(c % material) % rgb(2)
- b = pl % colors(c % material) % rgb(3)
- end if
- else if (pl % color_by == PLOT_COLOR_CELLS) then
- ! Assign color based on cell
- r = pl % colors(p % coord % cell) % rgb(1)
- g = pl % colors(p % coord % cell) % rgb(2)
- b = pl % colors(p % coord % cell) % rgb(3)
- else
- r = 0
- g = 0
- b = 0
- end if
- end if
+ ! get pixel color
+ call position_rgb(pl, rgb, id)
! Create a pixel at (x,y) with color (r,g,b)
- call set_pixel(img, x, y, r, g, b)
+ call set_pixel(img, x, y, rgb(1), rgb(2), rgb(3))
! Advance pixel in first direction
p % coord0 % xyz(in_i) = p % coord0 % xyz(in_i) + in_pixel
@@ -163,7 +178,7 @@ contains
subroutine output_ppm(pl, img)
- type(PlotSlice), pointer :: pl
+ type(ObjectPlot), pointer :: pl
type(Image), intent(in) :: img
integer :: i ! loop index for height
@@ -190,4 +205,81 @@ contains
end subroutine output_ppm
+!===============================================================================
+! CREATE_3D_DUMP outputs a binary file that can be input into silomesh for 3D
+! geometry visualization. It works the same way as create_ppm by dragging a
+! particle across the geometry for the specified number of voxels. The first
+! 3 int(4)'s in the binary are the number of x, y, and z voxels. The next 3
+! real(8)'s are the widths of the voxels in the x, y, and z directions. The next
+! 3 real(8)'s are the x, y, and z coordinates of the lower left point. Finally
+! the binary is filled with entries of four int(4)'s each. Each 'row' in the
+! binary contains four int(4)'s: 3 for x,y,z position and 1 for cell or material
+! id. For 1 million voxels this produces a file of approximately 15MB.
+!===============================================================================
+
+ subroutine create_3d_dump(pl)
+
+ type(ObjectPlot), pointer :: pl
+
+ integer :: x, y, z ! voxel location indices
+ integer :: rgb(3) ! colors (red, green, blue) from 0-255
+ integer :: id ! id of cell or material
+ real(8) :: vox(3) ! x, y, and z voxel widths
+ real(8) :: ll(3) ! lower left starting point for each sweep direction
+
+ ! compute voxel widths in each direction
+ vox = pl % width/dble(pl % pixels)
+
+ ! initial particle position
+ ll = pl % origin - pl % width / 2.0
+
+ ! allocate and initialize particle
+ allocate(p)
+ call initialize_particle()
+ p % coord0 % xyz = ll
+ p % coord0 % uvw = (/ 0.5, 0.5, 0.5 /)
+ p % coord0 % universe = BASE_UNIVERSE
+
+ ! Open binary plot file for writing
+ open(UNIT=UNIT_PLOT, FILE=pl % path_plot, STATUS='replace', &
+ ACCESS='stream')
+
+ ! write plot header info
+ write(UNIT_PLOT) pl % pixels, vox, ll
+
+ ! move to center of voxels
+ ll = ll + vox / 2.0
+
+ do x = 1, pl % pixels(1)
+ do y = 1, pl % pixels(2)
+ do z = 1, pl % pixels(3)
+
+ ! get voxel color
+ call position_rgb(pl, rgb, id)
+
+ ! write to plot file
+ write(UNIT_PLOT) id
+
+ ! advance particle in z direction
+ p % coord0 % xyz(3) = p % coord0 % xyz(3) + vox(3)
+
+ end do
+
+ ! advance particle in y direction
+ p % coord0 % xyz(2) = p % coord0 % xyz(2) + vox(2)
+ p % coord0 % xyz(3) = ll(3)
+
+ end do
+
+ ! advance particle in y direction
+ p % coord0 % xyz(1) = p % coord0 % xyz(1) + vox(1)
+ p % coord0 % xyz(2) = ll(2)
+ p % coord0 % xyz(3) = ll(3)
+
+ end do
+
+ close(UNIT_PLOT)
+
+ end subroutine create_3d_dump
+
end module plot
diff --git a/src/plot_header.F90 b/src/plot_header.F90
index 31ab15aa79..6956bfaf93 100644
--- a/src/plot_header.F90
+++ b/src/plot_header.F90
@@ -16,7 +16,7 @@ module plot_header
! PLOTSLICE holds plot information
!===============================================================================
- type PlotSlice
+ type ObjectPlot
integer :: id ! Unique ID
character(MAX_LINE_LEN) :: path_plot ! path for plot file
integer :: type ! Type
@@ -24,14 +24,14 @@ module plot_header
real(8) :: origin(3) ! xyz center of plot location
real(8) :: width(3) ! xyz widths of plot
integer :: basis ! direction of plot slice
- integer :: pixels(2) ! pixel width/height of plot slice
+ integer :: pixels(3) ! pixel width/height of plot slice
type(ObjectColor) :: not_found ! color for positions where no cell found
type(ObjectColor), allocatable :: colors(:) ! colors of cells/mats
- end type PlotSlice
+ end type ObjectPlot
- ! Plot type -- note that only slice plots are implemented currently
+ ! Plot type
integer, parameter :: PLOT_TYPE_SLICE = 1
- integer, parameter :: PLOT_TYPE_POINTS = 2
+ integer, parameter :: PLOT_TYPE_VOXEL = 2
! Plot basis plane
integer, parameter :: PLOT_BASIS_XY = 1
diff --git a/src/utils/plot_mesh_tally.py b/src/utils/plot_mesh_tally.py
new file mode 100755
index 0000000000..a759f7187c
--- /dev/null
+++ b/src/utils/plot_mesh_tally.py
@@ -0,0 +1,296 @@
+#!/usr/bin/env python
+
+'''Python script to plot tally data generated by OpenMC.'''
+
+import sys
+from statepoint import *
+
+# Color intensity dependent on individual score?
+
+from PyQt4.QtCore import *
+from PyQt4.QtGui import *
+import matplotlib.pyplot as plt
+from matplotlib.figure import Figure
+from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas
+from matplotlib.backends.backend_qt4agg import NavigationToolbar2QTAgg as NavigationToolbar
+import numpy as np
+
+class AppForm(QMainWindow):
+ def __init__(self, parent=None):
+ QMainWindow.__init__(self, parent)
+
+ # Read data from source or leakage fraction file
+ self.get_file_data()
+ self.main_frame = QWidget()
+ self.setCentralWidget(self.main_frame)
+
+ # Create the Figure, Canvas, and Axes
+ self.dpi = 100
+ self.fig = Figure((5.0, 15.0), dpi=self.dpi)
+ self.canvas = FigureCanvas(self.fig)
+ self.canvas.setParent(self.main_frame)
+ self.axes = self.fig.add_subplot(111)
+
+ # Create the navigation toolbar, tied to the canvas
+ self.mpl_toolbar = NavigationToolbar(self.canvas, self.main_frame)
+
+ # Grid layout at bottom
+ self.grid = QGridLayout()
+
+ # Overall layout
+ self.vbox = QVBoxLayout()
+ self.vbox.addWidget(self.canvas)
+ self.vbox.addWidget(self.mpl_toolbar)
+ self.vbox.addLayout(self.grid)
+ self.main_frame.setLayout(self.vbox)
+
+ # Tally selections
+ label_tally = QLabel("Tally:")
+ self.tally = QComboBox()
+ self.tally.addItems([(str(i + 1)) for i in range(self.n_tallies)])
+ self.connect(self.tally, SIGNAL('activated(int)'),
+ self._update)
+ self.connect(self.tally, SIGNAL('activated(int)'),
+ self.populate_boxes)
+ self.connect(self.tally, SIGNAL('activated(int)'),
+ self.on_draw)
+
+ # Planar basis
+ label_basis = QLabel("Basis:")
+ self.basis = QComboBox()
+ self.basis.addItems(['xy', 'yz', 'xz'])
+
+ # Update window when 'Basis' selection is changed
+ self.connect(self.basis, SIGNAL('activated(int)'),
+ self._update)
+ self.connect(self.basis, SIGNAL('activated(int)'),
+ self.populate_boxes)
+ self.connect(self.basis, SIGNAL('activated(int)'),
+ self.on_draw)
+
+ # Axial level within selected basis
+ label_axial_level = QLabel("Axial Level:")
+ self.axial_level = QComboBox()
+ self.connect(self.axial_level, SIGNAL('activated(int)'),
+ self.on_draw)
+
+ self.label_filters = QLabel("Filter options:")
+
+ # Labels for all possible filters
+ self.labels = {'cell': 'Cell: ', 'cellborn': 'Cell born: ',
+ 'surface': 'Surface: ', 'material': 'Material',
+ 'universe': 'Universe: ', 'energyin': 'Energy in: ',
+ 'energyout': 'Energy out: '}
+
+ # Empty reusable labels
+ self.qlabels = {}
+ for j in range(8):
+ self.nextLabel = QLabel
+ self.qlabels[j] = self.nextLabel
+
+ # Reusable comboboxes labelled with filter names
+ self.boxes = {}
+ for key in self.labels.keys():
+ self.nextBox = QComboBox()
+ self.connect(self.nextBox, SIGNAL('activated(int)'),
+ self.on_draw)
+ self.boxes[key] = self.nextBox
+
+ # Combobox to select among scores
+ self.score_label = QLabel("Score:")
+ self.scoreBox = QComboBox()
+ for item in self.tally_scores[0]:
+ self.scoreBox.addItems(str(item))
+ self.connect(self.scoreBox, SIGNAL('activated(int)'),
+ self.on_draw)
+
+ # Fill layout
+ self.grid.addWidget(label_tally, 0, 0)
+ self.grid.addWidget(self.tally, 0, 1)
+ self.grid.addWidget(label_basis, 1, 0)
+ self.grid.addWidget(self.basis, 1, 1)
+ self.grid.addWidget(label_axial_level, 2, 0)
+ self.grid.addWidget(self.axial_level, 2, 1)
+ self.grid.addWidget(self.label_filters, 3, 0)
+
+ self._update()
+ self.populate_boxes()
+ self.on_draw()
+
+ def get_file_data(self):
+ # Get data file name from "open file" browser
+ filename = QFileDialog.getOpenFileName(self, 'Select statepoint file', '.')
+
+ # Create StatePoint object and read in data
+ self.datafile = StatePoint(str(filename))
+ self.datafile.read_results()
+ self.datafile.generate_stdev()
+
+ self.setWindowTitle('Core Map Tool : ' + str(self.datafile.path))
+
+ # Set maximum colorbar value by maximum tally data value
+ self.maxvalue = self.datafile.tallies[0].results.max()
+
+ self.labelList = []
+
+ # Read mesh dimensions
+# for mesh in self.datafile.meshes:
+# self.nx, self.ny, self.nz = mesh.dimension
+
+ # Read filter types from statepoint file
+ self.n_tallies = len(self.datafile.tallies)
+ self.tally_list = []
+ for tally in self.datafile.tallies:
+ self.filter_types = []
+ for f in tally.filters:
+ self.filter_types.append(f)
+ self.tally_list.append(self.filter_types)
+
+ # Read score types from statepoint file
+ self.tally_scores = []
+ for tally in self.datafile.tallies:
+ self.score_types = []
+ for s in tally.scores:
+ self.score_types.append(s)
+ self.tally_scores.append(self.score_types)
+# print 'self.tally_scores = ', self.tally_scores
+
+ def on_draw(self):
+ """ Redraws the figure
+ """
+
+# print 'Calling on_draw...'
+ # Get selected basis, axial_level and stage
+ basis = self.basis.currentIndex() + 1
+ axial_level = self.axial_level.currentIndex() + 1
+
+ # Create spec_list
+ spec_list = []
+ for tally in self.datafile.tallies[self.tally.currentIndex()].filters.values():
+ if tally.type == 'mesh':
+ continue
+ index = self.boxes[tally.type].currentIndex()
+ spec_list.append((tally.type, index))
+
+ if self.basis.currentText() == '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(self.tally.currentIndex(), spec_list + [('mesh', (i, j, axial_level))], self.scoreBox.currentIndex())[0]
+
+ elif self.basis.currentText() == '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(self.tally.currentIndex(), spec_list + [('mesh', (axial_level, i, j))], self.scoreBox.currentIndex())[0]
+
+ 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(self.tally.currentIndex(), spec_list + [('mesh', (i, axial_level, j))], self.scoreBox.currentIndex())[0]
+
+# print spec_list
+
+ # Clear the figure
+ self.fig.clear()
+
+ # Make figure, set up color bar
+ self.axes = self.fig.add_subplot(111)
+ cax = self.axes.imshow(matrix, vmin=0.0, vmax=matrix.max(), interpolation="nearest")
+ self.fig.colorbar(cax)
+
+ self.axes.set_xticks([])
+ self.axes.set_yticks([])
+ self.axes.set_aspect('equal')
+
+ # Draw canvas
+ self.canvas.draw()
+
+ def _update(self):
+ '''Updates widget to display new relevant comboboxes and figure data
+ '''
+# print 'Calling _update...'
+
+ self.mesh = self.datafile.meshes[self.datafile.tallies[self.tally.currentIndex()].filters['mesh'].bins[0] - 1]
+
+ self.nx, self.ny, self.nz = self.mesh.dimension
+
+ # Clear axial level combobox
+ self.axial_level.clear()
+
+ # Repopulate axial level combobox based on current basis selection
+ if (self.basis.currentText() == 'xy'):
+ self.axial_level.addItems([str(i+1) for i in range(self.nz)])
+ elif (self.basis.currentText() == 'yz'):
+ self.axial_level.addItems([str(i+1) for i in range(self.nx)])
+ else:
+ self.axial_level.addItems([str(i+1) for i in range(self.ny)])
+
+ # Determine maximum value from current tally data set
+ self.maxvalue = self.datafile.tallies[self.tally.currentIndex()].results.max()
+# print self.maxvalue
+
+ # Clear and hide old filter labels
+ for item in self.labelList:
+ item.clear()
+
+ # Clear and hide old filter boxes
+ for j in self.labels:
+ self.boxes[j].clear()
+ self.boxes[j].setParent(None)
+
+ self.update()
+
+ def populate_boxes(self):
+# print 'Calling populate_boxes...'
+
+ n = 4
+ labels = {'cell': 'Cell : ',
+ 'cellborn': 'Cell born: ',
+ 'surface': 'Surface: ',
+ 'material': 'Material: ',
+ 'universe': 'Universe: '}
+
+ # For each filter in newly-selected tally, name a label and fill the
+ # relevant combobox with options
+ for element in self.tally_list[self.tally.currentIndex()]:
+ nextFilter = self.datafile.tallies[self.tally.currentIndex()].filters[element]
+ if element == 'mesh':
+ continue
+
+ label = QLabel(self.labels[element])
+ self.labelList.append(label)
+ combobox = self.boxes[element]
+ self.grid.addWidget(label, n, 0)
+ self.grid.addWidget(combobox, n, 1)
+ n += 1
+
+# print element
+ if element in ['cell', 'cellborn', 'surface', 'material', 'universe']:
+ combobox.addItems([str(i) for i in nextFilter.bins])
+# for i in nextFilter.bins:
+# print i
+
+ elif element == 'energyin' or element == 'energyout':
+ for i in range(nextFilter.length):
+ text = str(nextFilter.bins[i]) + ' to ' + str(nextFilter.bins[i+1])
+ combobox.addItem(text)
+
+ self.scoreBox.clear()
+ for item in self.tally_scores[self.tally.currentIndex()]:
+ self.scoreBox.addItem(str(item))
+ self.grid.addWidget(self.score_label, n, 0)
+ self.grid.addWidget(self.scoreBox, n, 1)
+
+
+
+def main():
+ app = QApplication(sys.argv)
+ form = AppForm()
+ form.show()
+ app.exec_()
+
+
+if __name__ == "__main__":
+ main()
diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py
index e91d474254..b6c23ea90c 100644
--- a/src/utils/statepoint.py
+++ b/src/utils/statepoint.py
@@ -555,30 +555,30 @@ class StatePoint(object):
def _get_int(self, n=1, path=None):
if self._hdf5:
- return self._f[path].value
+ return [int(v) for v in self._f[path].value]
else:
- return self._get_data(n, 'i', 4)
+ return [int(v) for v in self._get_data(n, 'i', 4)]
def _get_long(self, n=1, path=None):
if self._hdf5:
- return self._f[path].value
+ return [long(v) for v in self._f[path].value]
else:
- return self._get_data(n, 'q', 8)
+ return [long(v) for v in self._get_data(n, 'q', 8)]
def _get_float(self, n=1, path=None):
if self._hdf5:
- return self._f[path].value
+ return [float(v) for v in self._f[path].value]
else:
- return self._get_data(n, 'f', 4)
+ return [float(v) for v in self._get_data(n, 'f', 4)]
def _get_double(self, n=1, path=None):
if self._hdf5:
- return self._f[path].value
+ return [float(v) for v in self._f[path].value]
else:
- return self._get_data(n, 'd', 8)
+ return [float(v) for v in self._get_data(n, 'd', 8)]
def _get_string(self, n=1, path=None):
if self._hdf5:
- return self._f[path].value
+ return str(self._f[path].value)
else:
- return self._get_data(n, 's', 1)[0]
+ return str(self._get_data(n, 's', 1)[0])
diff --git a/src/utils/statepoint_silomesh.py b/src/utils/statepoint_3d.py
similarity index 75%
rename from src/utils/statepoint_silomesh.py
rename to src/utils/statepoint_3d.py
index f40f2ac1c9..13e5c10033 100755
--- a/src/utils/statepoint_silomesh.py
+++ b/src/utils/statepoint_3d.py
@@ -1,27 +1,36 @@
#!/usr/bin/env python
+from __future__ import division
+
+import sys
import itertools
import re
import warnings
-import silomesh # https://github.com/nhorelik/silomesh/
-
from statepoint import StatePoint
alphanum = re.compile(r"[\W_]+")
+err = False
+
################################################################################
def parse_options():
"""Process command line arguments"""
+
+
def tallies_callback(option, opt, value, parser):
"""Option parser function for list of tallies"""
+ global err
try:
setattr(parser.values, option.dest, [int(v) for v in value.split(',')])
- except: p.print_help()
+ except:
+ p.print_help()
+ err = True
def scores_callback(option, opt, value, parser):
"""Option parser function for list of scores"""
+ global err
try:
scores = {}
entries = value.split(',')
@@ -30,10 +39,13 @@ def parse_options():
if not tally in scores: scores[tally] = []
scores[tally].append(score)
setattr(parser.values, option.dest, scores)
- except: p.print_help()
+ except:
+ p.print_help()
+ err = True
def filters_callback(option, opt, value, parser):
"""Option parser function for list of filters"""
+ global err
try:
filters = {}
entries = value.split(',')
@@ -45,17 +57,17 @@ def parse_options():
filters[tally][filter_].append(bin)
setattr(parser.values, option.dest, filters)
except:
- raise
p.print_help()
+ err = True
from optparse import OptionParser
usage = r"""%prog [options]
-The default is to process all tallies and all scores into one silo file. Subsets
+The default is to process all tallies and all scores into one file. Subsets
can be chosen using the options. For example, to only process tallies 2 and 4
with all scores on tally 2 and only scores 1 and 3 on tally 4:
-%prog -t 2,4 -f 4.1,4.3
+%prog -t 2,4 -s 4.1,4.3
Likewise if you have additional filters on a tally you can specify a subset of
bins for each filter for that tally. For example to process all tallies and
@@ -72,34 +84,36 @@ You can list the available tallies, scores, and filters with the -l option:
help='List of tally indices to process, separated by commas.' \
' Default is to process all tallies.')
p.add_option('-s', '--scores', dest='scores', type='string', default=None,
- action='callback', callback=scores_callback,
- help='List of score indices to process, separated by commas, ' \
+ action='callback', callback=scores_callback,
+ help='List of score indices to process, separated by commas, ' \
'specified as {tallyid}.{scoreid}.' \
' Default is to process all scores in each tally.')
p.add_option('-f', '--filters', dest='filters', type='string', default=None,
- action='callback', callback=filters_callback,
- help='List of filter bins to process, separated by commas, ' \
+ action='callback', callback=filters_callback,
+ help='List of filter bins to process, separated by commas, ' \
'specified as {tallyid}.{filter}.{binid}. ' \
'Default is to process all filter combinaiton for each score.')
p.add_option('-l', '--list', dest='list', action='store_true',
help='List the tally and score indices available in the file.')
p.add_option('-o', '--output', action='store', dest='output',
- default='tally.silo', help='path to output SILO file.')
+ default='tally', help='path to output SILO file.')
p.add_option('-e', '--error', dest='valerr', default=False,
action='store_true',
help='Flag to extract errors instead of values.')
+ p.add_option('-v', '--vtk', action='store_true', dest='vtk',
+ default=False, help='Flag to convert to VTK instead of SILO.')
parsed = p.parse_args()
if not parsed[1]:
p.print_help()
- return parsed
+ return parsed, err
if parsed[0].valerr:
parsed[0].valerr = 1
else:
parsed[0].valerr = 0
- return parsed
+ return parsed, err
################################################################################
def main(file_, o):
@@ -113,8 +127,34 @@ def main(file_, o):
if o.list:
print_available(sp)
return
+
+ if o.vtk:
+ if not o.output[-4:] == ".vtm": o.output += ".vtm"
+ else:
+ if not o.output[-5:] == ".silo": o.output += ".silo"
+
+ if o.vtk:
+ try:
+ import vtk
+ except:
+ print 'The vtk python bindings do not appear to be installed properly.\n'+\
+ 'On Ubuntu: sudo apt-get install python-vtk\n'+\
+ 'See: http://www.vtk.org/'
+ return
+ else:
+ try:
+ import silomesh
+ except:
+ print 'The silomesh package does not appear to be installed properly.\n'+\
+ 'See: https://github.com/nhorelik/silomesh/'
+ return
- silomesh.init_silo(o.output)
+ if o.vtk:
+ blocks = vtk.vtkMultiBlockDataSet()
+ blocks.SetNumberOfBlocks(5)
+ block_idx = 0
+ else:
+ silomesh.init_silo(o.output)
# Tally loop #################################################################
for tally in sp.tallies:
@@ -128,8 +168,18 @@ def main(file_, o):
# extract filter options and mesh parameters for this tally
filtercombos = get_filter_combos(tally)
meshparms = get_mesh_parms(sp, tally)
- nx,ny,nz = meshparms[0], meshparms[1], meshparms[2]
- silomesh.init_mesh('Tally_{}'.format(tally.id), *meshparms)
+ nx,ny,nz = meshparms[:3]
+ ll = meshparms[3:6]
+ ur = meshparms[6:9]
+
+ if o.vtk:
+ ww = [(u-l)/n for u,l,n in zip(ur,ll,(nx,ny,nz))]
+ grid = grid = vtk.vtkImageData()
+ grid.SetDimensions(nx+1,ny+1,nz+1)
+ grid.SetOrigin(*ll)
+ grid.SetSpacing(*ww)
+ else:
+ silomesh.init_mesh('Tally_{}'.format(tally.id), *meshparms)
# Score loop ###############################################################
for sid,score in enumerate(tally.scores):
@@ -153,28 +203,58 @@ def main(file_, o):
# find and sanitize the variable name for this score
varname = get_sanitized_filterspec_name(tally, score, filterspec)
- silomesh.init_var(varname)
+ if o.vtk:
+ vtkdata = vtk.vtkDoubleArray()
+ vtkdata.SetName(varname)
+ dataforvtk = {}
+ else:
+ silomesh.init_var(varname)
- print "\t Score {}.{} {}:\t\t{}".format(tally.id, sid+1, score, varname)
+ lbl = "\t Score {}.{} {}:\t\t{}".format(tally.id, sid+1, score, varname)
# Mesh fill loop #######################################################
for x in range(1,nx+1):
+ sys.stdout.write(lbl+" {0}%\r".format(int(x/nx*100)))
+ sys.stdout.flush()
for y in range(1,ny+1):
for z in range(1,nz+1):
filterspec[0][1] = (x,y,z)
val = sp.get_value(tally.id-1, filterspec, sid)[o.valerr]
- silomesh.set_value(float(val), x, y, z)
-
+ if o.vtk:
+ # vtk cells go z, y, x, so we store it now and enter it later
+ i = (z-1)*nx*ny + (y-1)*nx + x-1
+ dataforvtk[i] = float(val)
+ else:
+ silomesh.set_value(float(val), x, y, z)
+
# end mesh fill loop
- silomesh.finalize_var()
+ print
+ if o.vtk:
+ for i in range(nx*ny*nz):
+ vtkdata.InsertNextValue(dataforvtk[i])
+ grid.GetCellData().AddArray(vtkdata)
+ del vtkdata
+
+ else:
+ silomesh.finalize_var()
# end filter loop
# end score loop
- silomesh.finalize_mesh()
+ if o.vtk:
+ blocks.SetBlock(block_idx, grid)
+ block_idx += 1
+ else:
+ silomesh.finalize_mesh()
# end tally loop
- silomesh.finalize_silo()
+ if o.vtk:
+ writer = vtk.vtkXMLMultiBlockDataWriter()
+ writer.SetFileName(o.output)
+ writer.SetInput(blocks)
+ writer.Write()
+ else:
+ silomesh.finalize_silo()
################################################################################
def get_sanitized_filterspec_name(tally, score, filterspec):
@@ -237,7 +317,7 @@ def print_available(sp):
mesh = ""
if not 'mesh' in tally.filters: mesh = "(no mesh)"
print "\tTally {} {}".format(tally.id, mesh)
- scores = ["{}.{}: {}".format(tally.id, sid+1, score)
+ scores = ["{}.{}: {}".format(tally.id, sid, score)
for sid, score in enumerate(tally.scores)]
for score in scores:
print "\t\tScore {}".format(score)
@@ -250,8 +330,8 @@ def print_available(sp):
def validate_options(sp,o):
"""Validates specified tally/score options for the current statepoint"""
+ available_tallies = [t.id for t in sp.tallies]
if o.tallies:
- available_tallies = [t.id for t in sp.tallies]
for otally in o.tallies:
if not otally in available_tallies:
warnings.warn('Tally {} not in statepoint file'.format(otally))
@@ -308,6 +388,6 @@ warnings.formatwarning = formatwarning
################################################################################
if __name__ == '__main__':
- (options, args) = parse_options()
- if args:
+ (options, args), err = parse_options()
+ if args and not err:
main(args[0],options)
diff --git a/src/utils/voxel.py b/src/utils/voxel.py
new file mode 100755
index 0000000000..c0d7e78d3d
--- /dev/null
+++ b/src/utils/voxel.py
@@ -0,0 +1,120 @@
+#!/usr/bin/env python
+
+from __future__ import division
+
+import struct
+import sys
+
+################################################################################
+def parse_options():
+ """Process command line arguments"""
+
+ from optparse import OptionParser
+ usage = r"""%prog [options] """
+ p = OptionParser(usage=usage)
+ p.add_option('-o', '--output', action='store', dest='output',
+ default='plot', help='Path to output SILO or VTK file.')
+ p.add_option('-v', '--vtk', action='store_true', dest='vtk',
+ default=False, help='Flag to convert to VTK instead of SILO.')
+ parsed = p.parse_args()
+ if not parsed[1]:
+ p.print_help()
+ return parsed
+ return parsed
+
+################################################################################
+def main(file_, o):
+
+ print file_
+ fh = open(file_,'rb')
+ header = get_header(fh)
+ meshparms = header['dimension'] + header['lower_left'] + header['upper_right']
+ nx,ny,nz = meshparms[0], meshparms[1], meshparms[2]
+ ll = header['lower_left']
+
+ if o.vtk:
+ try:
+ import vtk
+ except:
+ print 'The vtk python bindings do not appear to be installed properly.\n'+\
+ 'On Ubuntu: sudo apt-get install python-vtk\n'+\
+ 'See: http://www.vtk.org/'
+ return
+
+ origin = [(l+w*n/2.) for n,l,w in zip((nx,ny,nz),ll,header['width'])]
+
+ grid = vtk.vtkImageData()
+ grid.SetDimensions(nx+1,ny+1,nz+1)
+ grid.SetOrigin(*ll)
+ grid.SetSpacing(*header['width'])
+
+ data = vtk.vtkDoubleArray()
+ data.SetName("id")
+ data.SetNumberOfTuples(nx*ny*nz)
+ for x in range(nx):
+ sys.stdout.write(" {0}%\r".format(int(x/nx*100)))
+ sys.stdout.flush()
+ for y in range(ny):
+ for z in range(nz):
+ i = z*nx*ny + y*nx + x
+ id_ = get_int(fh)[0]
+ data.SetValue(i, id_)
+ grid.GetCellData().AddArray(data)
+
+ writer = vtk.vtkXMLImageDataWriter()
+ writer.SetInput(grid)
+ if not o.output[-4:] == ".vti": o.output += ".vti"
+ writer.SetFileName(o.output)
+ writer.Write()
+
+ else:
+
+ try:
+ import silomesh
+ except:
+ print 'The silomesh package does not appear to be installed properly.\n'+\
+ 'See: https://github.com/nhorelik/silomesh/'
+ return
+ if not o.output[-5:] == ".silo": o.output += ".silo"
+ silomesh.init_silo(o.output)
+ silomesh.init_mesh('plot', *meshparms)
+ silomesh.init_var("id")
+ for x in range(1,nx+1):
+ sys.stdout.write(" {0}%\r".format(int(x/nx*100)))
+ sys.stdout.flush()
+ for y in range(1,ny+1):
+ for z in range(1,nz+1):
+ id_ = get_int(fh)[0]
+ silomesh.set_value(float(id_), x, y, z)
+ print
+ silomesh.finalize_var()
+ silomesh.finalize_mesh()
+ silomesh.finalize_silo()
+
+################################################################################
+def get_header(file_):
+ nx,ny,nz = get_int(file_, 3)
+ wx,wy,wz = get_double(file_, 3)
+ lx,ly,lz = get_double(file_, 3)
+ header = {'dimension':[nx,ny,nz], 'width':[wx,wy,wz], 'lower_left':[lx,ly,lz],
+ 'upper_right': [lx+wx*nx,ly+wy*ny,lz+wz*nz]}
+ return header
+
+################################################################################
+def get_data(file_, n, typeCode, size):
+ return list(struct.unpack('={0}{1}'.format(n,typeCode),
+ file_.read(n*size)))
+
+################################################################################
+def get_int(file_, n=1, path=None):
+ return get_data(file_, n, 'i', 4)
+
+################################################################################
+def get_double(file_, n=1, path=None):
+ return get_data(file_, n, 'd', 8)
+
+################################################################################
+if __name__ == '__main__':
+ (options, args) = parse_options()
+ if args:
+ main(args[0],options)