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

Merge branch 'master' into HEAD

Browse source
This commit is contained in:
Paul Romano 2013-04-22 18:02:59 -04:00
commit 0ea4ba1579
24 changed files with 1295 additions and 191 deletions
Download patch file Download diff file Expand all files Collapse all files

BIN
docs/img/3dcore.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 278 KiB

BIN
docs/img/3dgeomplot.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 80 KiB

BIN
docs/img/atr.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 460 KiB

BIN
docs/img/fluxplot.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 26 KiB

BIN
docs/img/plotmeshtally.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 71 KiB

1
docs/source/devguide/index.rst
View file

@ -17,3 +17,4 @@ as debugging.
workflow
xml-fortran
statepoint
voxel

54
docs/source/devguide/voxel.rst Normal file
View file

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

4
docs/source/usersguide/beginners.rst
View file

@ -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
<http://www.ee.surrey.ac.uk/Teaching/Unix/>`_ will help you get acquainted with
commonly-used commands.
commonly-used commands. It is also helpful to be familiar with `Python
<http://www.python.org/>`_, 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

1
docs/source/usersguide/index.rst
View file

@ -14,4 +14,5 @@ essential aspects of using OpenMC to perform neutronic simulations.
beginners
install
input
processing
troubleshoot

114
docs/source/usersguide/input.rst
View file

@ -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 ``<plots>`` and any number output figures can
be defined with ``<plot>`` sub-elements.
element of the plots.xml is simply ``<plots>`` and any number output plots can
be defined with ``<plot>`` 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.
``<plot>`` 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"
``<plot>`` elements of ``type`` "slice" also contain the following attributes or
sub-elements:
``<plot>`` 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
``<plot>`` 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 ``<plot>`` 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 ``<plot>`` 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

360
docs/source/usersguide/processing.rst Normal file
View file

@ -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 <http://www.numpy.org/>`_
* [1]_ `Scipy <http://www.scipy.org/>`_
* [2]_ `h5py <http://code.google.com/p/h5py/>`_
* [3]_ `Matplotlib <http://matplotlib.org/>`_
* [3]_ `Silomesh <https://github.com/nhorelik/silomesh>`_
* [3]_ `VTK <http://www.vtk.org/>`_
* [4]_ `PyQt <http://www.riverbankcomputing.com/software/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 <http://www.gimp.org/>`_, `IrfanView
<http://www.irfanview.com/>`_, 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 <http://www.imagemagick.org/script/convert.php>`_ 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
<openmc_root>/src/utils/voxel.py myplot.voxel -o output.silo
and VTK file formats:
.. code-block:: sh
<openmc_root>/src/utils/myplot.voxel --vtk -o output.vti
To use this utility you need either
* `Silomesh <https://github.com/nhorelik/silomesh>`_
or
* `VTK <http://www.vtk.org/>`_ 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 <http://ipython.org/>`_). 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 <http://www.riverbankcomputing.com/software/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
<openmc_root>/src/utils/statepoint_3d.py <statepoint_file> -o output.silo
<openmc_root>/src/utils/statepoint_3d.py <statepoint_file> --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
<openmc_root>/src/utils/statepoint_3d.py <statepoint_file> --tallies 2,4 --scores 4.1,4.3 -o output.silo
<openmc_root>/src/utils/statepoint_3d.py <statepoint_file> --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
<http://docs.scipy.org/doc/scipy/reference/tutorial/io.html>`_ 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.

13
src/constants.F90
View file

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

11
src/geometry.F90
View file

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

4
src/global.F90
View file

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

1
src/hdf5_interface.F90
View file

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

132
src/input_xml.F90
View file

@ -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 = "<pixels> 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 = "<pixels> 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 = "<pixels> 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 = "<width> 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 = "<width> 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" // &

31
src/output.F90
View file

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

2
src/physics.F90
View file

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

176
src/plot.F90
View file

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

10
src/plot_header.F90
View file

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

296
src/utils/plot_mesh_tally.py Executable file
View file

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

20
src/utils/statepoint.py
View file

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

136
src/utils/statepoint_silomesh.py → src/utils/statepoint_3d.py
View file

@ -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] <statepoint_file>
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 <statepoint_file>
%prog -t 2,4 -s 4.1,4.3 <statepoint_file>
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)

120
src/utils/voxel.py Executable file
View file

@ -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] <voxel_file>"""
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)