Merge branch 'master' into raster3d

This commit is contained in:
nhorelik 2013-04-13 22:15:54 -04:00
commit 3af332318c
15 changed files with 663 additions and 269 deletions

View file

@ -332,15 +332,14 @@ source.o: output.o
source.o: particle_header.o
source.o: physics.o
source.o: random_lcg.o
source.o: state_point.o
source.o: string.o
state_point.o: error.o
state_point.o: global.o
state_point.o: math.o
state_point.o: output.o
state_point.o: source.o
state_point.o: string.o
state_point.o: tally.o
state_point.o: tally_header.o
string.o: constants.o

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)

View file

@ -9,6 +9,7 @@ module source
use particle_header, only: deallocate_coord
use physics, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn, set_particle_seed
use state_point, only: read_source_binary
use string, only: to_str
#ifdef MPI
@ -229,175 +230,4 @@ contains
end subroutine initialize_particle
!===============================================================================
! WRITE_SOURCE writes out the final source distribution to a binary file that
! can be used as a starting source in a new simulation
!===============================================================================
subroutine write_source_binary()
#ifdef MPI
integer :: fh ! file handle
integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file)
! ==========================================================================
! PARALLEL I/O USING MPI-2 ROUTINES
! Open binary source file for reading
call MPI_FILE_OPEN(MPI_COMM_WORLD, path_source, MPI_MODE_CREATE + &
MPI_MODE_WRONLY, MPI_INFO_NULL, fh, mpi_err)
if (master) then
offset = 0
call MPI_FILE_WRITE_AT(fh, offset, n_particles, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpi_err)
end if
! Set proper offset for source data on this processor
offset = 8*(1 + rank*maxwork*8)
! Write all source sites
call MPI_FILE_WRITE_AT(fh, offset, source_bank(1), work, MPI_BANK, &
MPI_STATUS_IGNORE, mpi_err)
! Close binary source file
call MPI_FILE_CLOSE(fh, mpi_err)
#else
! ==========================================================================
! SERIAL I/O USING FORTRAN INTRINSIC ROUTINES
! Open binary source file for writing
open(UNIT=UNIT_SOURCE, FILE=path_source, STATUS='replace', &
ACCESS='stream')
! Write the number of particles
write(UNIT=UNIT_SOURCE) n_particles
! Write information from the source bank
write(UNIT=UNIT_SOURCE) source_bank(1:work)
! Close binary source file
close(UNIT=UNIT_SOURCE)
#endif
end subroutine write_source_binary
!===============================================================================
! READ_SOURCE_BINARY reads a source distribution from a source.binary file and
! initializes the source bank
!===============================================================================
subroutine read_source_binary()
integer :: i ! loop over repeating sites
integer(8) :: n_sites ! number of sites in binary file
integer :: n_repeat ! number of times to repeat a site
#ifdef MPI
integer :: fh ! file handle
integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file)
integer :: n_read ! number of sites to read on a single process
#endif
#ifdef MPI
! ==========================================================================
! PARALLEL I/O USING MPI-2 ROUTINES
! Open binary source file for reading
call MPI_FILE_OPEN(MPI_COMM_WORLD, path_source, MPI_MODE_RDONLY, &
MPI_INFO_NULL, fh, mpi_err)
! Read number of source sites in file
offset = 0
call MPI_FILE_READ_AT(fh, offset, n_sites, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpi_err)
if (n_particles > n_sites) then
! Determine number of sites to read and offset
if (rank <= mod(n_sites,int(n_procs,8)) - 1) then
n_read = int(n_sites/n_procs) + 1
offset = 8*(1 + rank*n_read*8)
else
n_read = int(n_sites/n_procs)
offset = 8*(1 + (rank*n_read + mod(n_sites,int(n_procs,8)))*8)
end if
! Read source sites
call MPI_FILE_READ_AT(fh, offset, source_bank(1), n_read, MPI_BANK, &
MPI_STATUS_IGNORE, mpi_err)
! Let's say we have 30 sites and we need to fill in 200. This do loop
! will fill in sites 31 - 180.
n_repeat = int(work / n_read)
do i = 1, n_repeat - 1
source_bank(i*n_read + 1:(i+1)*n_read) = &
source_bank((i-1)*n_read + 1:i*n_read)
end do
! This final statement would fill sites 181 - 200 in the above example.
if (mod(work, int(n_repeat*n_read,8)) > 0) then
source_bank(n_repeat*n_read + 1:work) = &
source_bank(1:work - n_repeat * n_read)
end if
else
! Set proper offset for source data on this processor
offset = 8*(1 + rank*maxwork*8)
! Read all source sites
call MPI_FILE_READ_AT(fh, offset, source_bank(1), work, MPI_BANK, &
MPI_STATUS_IGNORE, mpi_err)
! Close binary source file
call MPI_FILE_CLOSE(fh, mpi_err)
end if
#else
! ==========================================================================
! SERIAL I/O USING FORTRAN INTRINSIC ROUTINES
! Open binary source file for reading
open(UNIT=UNIT_SOURCE, FILE=path_source, STATUS='old', &
ACCESS='stream')
! Read number of source sites in file
read(UNIT=UNIT_SOURCE) n_sites
if (n_particles > n_sites) then
! The size of the source file is smaller than the number of particles we
! need. Thus, read all sites and then duplicate sites as necessary.
read(UNIT=UNIT_SOURCE) source_bank(1:n_sites)
! Let's say we have 300 sites and we need to fill in 1000. This do loop
! will fill in sites 301 - 900.
n_repeat = int(n_particles / n_sites)
do i = 1, n_repeat - 1
source_bank(i*n_sites + 1:(i+1)*n_sites) = &
source_bank((i-1)*n_sites + 1:i*n_sites)
end do
! This final statement would fill sites 901 - 1000 in the above example.
source_bank(n_repeat*n_sites + 1:n_particles) = &
source_bank(1:n_particles - n_repeat * n_sites)
else
! The size of the source file is bigger than or equal to the number of
! particles we need for one generation. Thus, we can just read as many
! sites as we need.
read(UNIT=UNIT_SOURCE) source_bank(1:n_particles)
end if
! Close binary source file
close(UNIT=UNIT_SOURCE)
#endif
end subroutine read_source_binary
end module source

View file

@ -22,10 +22,8 @@ module state_point
use global
use math, only: t_percentile
use output, only: write_message, print_batch_keff, time_stamp
use source, only: write_source_binary
use string, only: to_str
use tally_header, only: TallyObject
use tally, only: setup_active_usertallies
#ifdef MPI
use mpi
@ -1113,4 +1111,175 @@ contains
end subroutine replay_batch_history
!===============================================================================
! WRITE_SOURCE writes out the final source distribution to a binary file that
! can be used as a starting source in a new simulation
!===============================================================================
subroutine write_source_binary()
#ifdef MPI
integer :: fh ! file handle
integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file)
! ==========================================================================
! PARALLEL I/O USING MPI-2 ROUTINES
! Open binary source file for reading
call MPI_FILE_OPEN(MPI_COMM_WORLD, path_source, MPI_MODE_CREATE + &
MPI_MODE_WRONLY, MPI_INFO_NULL, fh, mpi_err)
if (master) then
offset = 0
call MPI_FILE_WRITE_AT(fh, offset, n_particles, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpi_err)
end if
! Set proper offset for source data on this processor
offset = 8*(1 + rank*maxwork*8)
! Write all source sites
call MPI_FILE_WRITE_AT(fh, offset, source_bank(1), work, MPI_BANK, &
MPI_STATUS_IGNORE, mpi_err)
! Close binary source file
call MPI_FILE_CLOSE(fh, mpi_err)
#else
! ==========================================================================
! SERIAL I/O USING FORTRAN INTRINSIC ROUTINES
! Open binary source file for writing
open(UNIT=UNIT_SOURCE, FILE=path_source, STATUS='replace', &
ACCESS='stream')
! Write the number of particles
write(UNIT=UNIT_SOURCE) n_particles
! Write information from the source bank
write(UNIT=UNIT_SOURCE) source_bank(1:work)
! Close binary source file
close(UNIT=UNIT_SOURCE)
#endif
end subroutine write_source_binary
!===============================================================================
! READ_SOURCE_BINARY reads a source distribution from a source.binary file and
! initializes the source bank
!===============================================================================
subroutine read_source_binary()
integer :: i ! loop over repeating sites
integer(8) :: n_sites ! number of sites in binary file
integer :: n_repeat ! number of times to repeat a site
#ifdef MPI
integer :: fh ! file handle
integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file)
integer :: n_read ! number of sites to read on a single process
#endif
#ifdef MPI
! ==========================================================================
! PARALLEL I/O USING MPI-2 ROUTINES
! Open binary source file for reading
call MPI_FILE_OPEN(MPI_COMM_WORLD, path_source, MPI_MODE_RDONLY, &
MPI_INFO_NULL, fh, mpi_err)
! Read number of source sites in file
offset = 0
call MPI_FILE_READ_AT(fh, offset, n_sites, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpi_err)
if (n_particles > n_sites) then
! Determine number of sites to read and offset
if (rank <= mod(n_sites,int(n_procs,8)) - 1) then
n_read = int(n_sites/n_procs) + 1
offset = 8*(1 + rank*n_read*8)
else
n_read = int(n_sites/n_procs)
offset = 8*(1 + (rank*n_read + mod(n_sites,int(n_procs,8)))*8)
end if
! Read source sites
call MPI_FILE_READ_AT(fh, offset, source_bank(1), n_read, MPI_BANK, &
MPI_STATUS_IGNORE, mpi_err)
! Let's say we have 30 sites and we need to fill in 200. This do loop
! will fill in sites 31 - 180.
n_repeat = int(work / n_read)
do i = 1, n_repeat - 1
source_bank(i*n_read + 1:(i+1)*n_read) = &
source_bank((i-1)*n_read + 1:i*n_read)
end do
! This final statement would fill sites 181 - 200 in the above example.
if (mod(work, int(n_repeat*n_read,8)) > 0) then
source_bank(n_repeat*n_read + 1:work) = &
source_bank(1:work - n_repeat * n_read)
end if
else
! Set proper offset for source data on this processor
offset = 8*(1 + rank*maxwork*8)
! Read all source sites
call MPI_FILE_READ_AT(fh, offset, source_bank(1), work, MPI_BANK, &
MPI_STATUS_IGNORE, mpi_err)
! Close binary source file
call MPI_FILE_CLOSE(fh, mpi_err)
end if
#else
! ==========================================================================
! SERIAL I/O USING FORTRAN INTRINSIC ROUTINES
! Open binary source file for reading
open(UNIT=UNIT_SOURCE, FILE=path_source, STATUS='old', &
ACCESS='stream')
! Read number of source sites in file
read(UNIT=UNIT_SOURCE) n_sites
if (n_particles > n_sites) then
! The size of the source file is smaller than the number of particles we
! need. Thus, read all sites and then duplicate sites as necessary.
read(UNIT=UNIT_SOURCE) source_bank(1:n_sites)
! Let's say we have 300 sites and we need to fill in 1000. This do loop
! will fill in sites 301 - 900.
n_repeat = int(n_particles / n_sites)
do i = 1, n_repeat - 1
source_bank(i*n_sites + 1:(i+1)*n_sites) = &
source_bank((i-1)*n_sites + 1:i*n_sites)
end do
! This final statement would fill sites 901 - 1000 in the above example.
source_bank(n_repeat*n_sites + 1:n_particles) = &
source_bank(1:n_particles - n_repeat * n_sites)
else
! The size of the source file is bigger than or equal to the number of
! particles we need for one generation. Thus, we can just read as many
! sites as we need.
read(UNIT=UNIT_SOURCE) source_bank(1:n_particles)
end if
! Close binary source file
close(UNIT=UNIT_SOURCE)
#endif
end subroutine read_source_binary
end module state_point

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