Add new plotting tools and add performance improvements (91cd8451) · Commits · Phathanapirom, Birdy / SubKit

SubKit/process/Hdf5Tools.py

+105 −56

Original line number	Diff line number	Diff line
		@@ -36,6 +36,50 @@ class Hdf5Interface(object):
		for pin in range(me.numPin):
		me.pinNumAxial[pin+1] = len(me.h5['CORE/mesh/PIN_{:05d}/axial_nodes'.format(pin+1)])

		# Create a mask for pin surface datasets, which can be used to extract relevant data
		# True means the data is irrelevant
		# For now we assum that dnbr will always be in the file (if there are pins in the model)
		# If we ever make it optional, will need to do a search for a surf dataset
		if 'STATE_0001/pin_surf_dnbr' in me.h5:
		me.surfMask = np.ones(me.h5['STATE_0001/pin_surf_dnbr'].shape, dtype=bool)
		for pin in range(1, me.numPin+1):
		nlev = me.getPinNumLev(pin)
		nsec = me.getPinNumSector(pin)
		me.surfMask[0:nlev, 0:nsec, pin-1] = False

		me.numCh = 1
		while True:
		if '/CORE/mesh/CHAN_{:05d}'.format(me.numCh) not in me.h5:
		me.numCh = me.numCh-1
		break
		else:
		me.numCh = me.numCh+1

		# Create a mask for channel datasets
		# Need one for scalar data and one for momentum data
		me.chanMomMask = np.ones(me.h5['STATE_0001/chan_temp'].shape, dtype=bool)
		me.chanScalarMask = np.ones(me.h5['STATE_0001/chan_pressure'].shape, dtype=bool)
		for chan in range(1, me.numCh+1):
		nlev = len(me.h5['CORE/mesh/CHAN_{:05d}/scalar_bounds'.format(chan)])
		me.chanMomMask[0:nlev, chan-1] = False
		me.chanScalarMask[0:nlev+1, chan-1] = False

		# Create a mask for pin internal data
		if 'STATE_0001/pin_temp' in me.h5:
		me.pinInternalMask = np.ones(me.h5['STATE_0001/pin_temp'].shape, dtype=bool)
		for pin in range(1, me.numPin+1):
		nlev = me.getPinNumLev(pin)
		nsec = me.getPinNumSector(pin)
		nrad = me.getPinNumRadial(pin)
		me.pinInternalMask[0:nlev, 0:nsec, 0:nrad, pin-1] = False

		# Create a mask for pin axial data
		if 'STATE_0001/pin_axial_centerline_temp' in me.h5:
		me.pinAxialMask = np.ones(me.h5['STATE_0001/pin_axial_centerline_temp'].shape, dtype=bool)
		for pin in range(1, me.numPin+1):
		nlev = me.getPinNumLev(pin)
		me.pinAxialMask[0:nlev, pin-1] = False

		# Location the axial locations of the section boundaries
		me.sectionBounds = [0.0]
		for section in list(range(1, me.getNumSections()+1)):
		@@ -50,13 +94,6 @@ class Hdf5Interface(object):
		if 'mesh_type' in dataset.attrs:
		me.chanDatasetMeshTypes[datasetName] = dataset.attrs['mesh_type'][0]

		me.numCh = 1
		while True:
		if '/CORE/mesh/CHAN_{:05d}'.format(me.numCh) not in me.h5:
		me.numCh = me.numCh-1
		break
		else:
		me.numCh = me.numCh+1

		def getNumPin(me):
		""" Returns the total number of pins in the model """
		@@ -164,6 +201,16 @@ class Hdf5Interface(object):
		me._pinIndexValid(pin)
		return me.pinNumAxial[pin]

		def getChanArea(me, ch):
		""" Returns the nominal channel flow area [m**2]

		Args:
		ch (int) : Channel index (1-based)

		"""
		me._chanIndexValid(ch)
		return me.h5['CORE/chan_area'][ch-1].value

		def getPinThetas(me, pin, level):
		""" Returns the azimuthal fractions for the passed pin and level

		@@ -342,6 +389,23 @@ class Hdf5Interface(object):
		me._sumChanDataset('chan_mdot_drp', state, chans, _level)
		return mdotTot/me.getCrossSectionalArea(section)

		def getChanDatasetMeshType(me, dataset):
		""" Returns the mesh type this channel dataset uses

		Args:
		dataset (str) : Dataset name

		Returns:
		meshType (str): Will be "momentum" or "scalar"
		"""
		me._datasetValid(dataset)
		assert('mesh_type' in me.h5['STATE_0001/{:s}'.format(dataset)].attrs)
		meshType = me.h5['STATE_0001/{:s}'.format(dataset)].attrs['mesh_type'][0]
		if meshType=='scalar_bounds':
		return "momentum"
		else:
		return "scalar"

		def getChanDataset(me, dataset, state, chan):
		""" Retrieve dataset for a channel from the h5 file

		@@ -450,20 +514,13 @@ class Hdf5Interface(object):
		pins = [pin]
		maxVal = None
		for state_ in states:
		for pin in pins:
		data = me.getPinAxialDataset(dataset, state_, pin)
		maxData = max(data)
		if maxVal is None:
		maxVal = maxData
		ma = np.ma.array(me.h5['STATE_{:04d}/{:s}'.format(state_, dataset)], mask=me.pinAxialMask)[:, np.array(pins)-1]
		maxval = ma.max()
		if maxVal is None or maxval>maxVal:
		maxVal = maxval
		maxState = state_
		maxPin = pin
		maxLevel = np.argmax(data)
		elif maxData>maxVal:
		maxVal = maxData
		maxState = state_
		maxPin = pin
		maxLevel = np.argmax(data)
		return maxState, maxPin, maxLevel, maxVal
		(maxLevel, maxPin) = np.unravel_index(ma.argmax(), ma.shape)
		return maxState, maxPin+1, maxLevel, maxVal

		def _getPinSurfaceBounding(me, dataset, state=None, operation='max'):
		""" Use to find the largest or smallest value and location in a pin surface dataset
		@@ -490,33 +547,22 @@ class Hdf5Interface(object):
		else:
		states = [state]
		assert(operation=='max' or operation=='min')
		if operation=='max':
		op=np.max
		loc=np.argmax
		def pastCurrentBound(currentBound, newVal):
		if currentBound>newVal:
		return True
		else:
		return False
		else:
		op=np.min
		loc=np.argmin
		def pastCurrentBound(currentBound, newVal):
		if newVal<currentBound:
		return True
		else:
		return False
		globalBound = None
		for state_ in states:
		for pin in range(1, me.getNumPin()+1):
		data = me.getPinSurfaceDataset(dataset, state_, pin)
		boundVal = op(data)
		if globalBound is None or pastCurrentBound(globalBound, boundVal):
		globalBound = boundVal
		ma = np.ma.array(me.h5['STATE_{:04d}/{:s}'.format(state_, dataset)], mask=me.surfMask)
		if operation=='max':
		maxval = ma.max()
		if globalBound is None or maxval>globalBound:
		globalBound = maxval
		boundState = state_
		(boundLevel, boundSector, boundPin) = np.unravel_index(ma.argmax(), ma.shape)
		else:
		minval = ma.min()
		if globalBound is None or minval<globalBound:
		globalBound = minval
		boundState = state_
		boundPin = pin
		(boundLevel, boundSector) = np.unravel_index(loc(data, axis=None), data.shape)
		return boundState, boundPin, boundLevel, boundSector, globalBound
		(boundLevel, boundSector, boundPin) = np.unravel_index(ma.argmin(), ma.shape)
		return boundState, boundPin+1, boundLevel, boundSector, globalBound

		def getPinWithMaxCenterlineTemp(me, state=None, pin=None):
		""" Finds the pin with the maximum centerline temperature.
		@@ -829,17 +875,20 @@ class Hdf5Interface(object):
		else:
		searchStates = list(range(1, me.getNumState()+1))

		if me.getChanDatasetMeshType(dataset)=="scalar":
		mask = me.chanScalarMask
		else:
		mask = me.chanMomMask

		for st in searchStates:
		for chan in range(1, me.getNumChan()+1):
		data = me.getChanDataset(dataset, st, chan)
		maxLocal = max(data)
		if maxVal is None or maxLocal > maxVal:
		ma = np.ma.array(me.h5['STATE_{:04d}/{:s}'.format(st, dataset)], mask=mask)
		maxval = ma.max()
		if maxVal is None or maxval>maxVal:
		maxVal = maxval
		maxState = st
		maxVal = maxLocal
		maxChan = chan
		maxLevel = np.argmax(data)
		(maxLevel, maxChan) = np.unravel_index(ma.argmax(), ma.shape)

		return maxState, maxChan, maxLevel, maxVal
		return maxState, maxChan+1, maxLevel, maxVal


		def getChanHeight(me, chan):

SubKit/process/PlotTools.py

+46 −15

Original line number	Diff line number	Diff line
		@@ -59,7 +59,7 @@ class PlotTools(object):
		radialMesh, radialTemp = me.h5obj.getRadialTemp(state_, pin_, level_, includeCenterline=True)

		# Make all indices 1-based for I/O
		location = "STATE_{:04d}_PIN_{:05d}_level_{:03d}".format(state_, pin_, level_+1)
		location = "S{:d}_P{:d}_L{:d}".format(state_, pin_, level_+1)

		if not figname:
		figname_ = location+"_pin_temp.png"
		@@ -105,7 +105,7 @@ class PlotTools(object):
		theta_ = theta

		# Make all indices 1-based for I/O
		location = "STATE_{:04d}_PIN_{:05d}_theta_{:02d}".format(state_, pin_, theta_+1)
		location = "S{:d}_P{:d}_T{:d}".format(state_, pin_, theta_+1)

		if not figname:
		figname_ = location+"_pin_dnbr.png"
		@@ -126,7 +126,7 @@ class PlotTools(object):
		plt.close(fig)

		def _chName(me, state, ch):
		return "STATE_{:04d}_CHAN_{:05d}".format(state, ch)
		return "S{:04d}_C{:05d}".format(state, ch)

		def _plotAxial(me, x, y, xname, figname, labels, title=None, yname="Axial [m]"):
		""" Makes an axial plot
		@@ -166,15 +166,16 @@ class PlotTools(object):
		plt.savefig(figname)
		plt.close(fig)

		def plotChAxialVoid(me, state=None, ch=None, figname=None):
		""" Makes an axial plot of void for a channel
		def _chanPlotMaxWrapper(me, dataset, plotYLabel, state=None, ch=None, figname=None):
		""" A wrapper for plotting channel data

		Args:
		state (int): State number to plot (1-based). If not present, selects state with maximum.
		ch (int): Chan number to plot (1-based). If not present, selects chan with maximum.
		figname (str): Name of the figure.
		Will default to plotting the max if state/channels are not provided.
		Also provides logic for handling single vs. multiple channels.

		Accepts "massflux_tot" as a dataset which does not appear in the file but can be calculated
		"""
		stateMax, chMax, levMax, voidMax = me.h5obj.getMaxChanData('chan_void_vap', state)

		stateMax, chMax, levMax, maxVal = me.h5obj.getMaxChanData(dataset, state)
		if state is None:
		state_ = stateMax
		else:
		@@ -184,20 +185,50 @@ class PlotTools(object):
		else:
		assert(isinstance(ch, int) or isinstance(ch, list))
		if isinstance(ch, int):
		ch_ = [list(set(ch))]
		ch_ = [ch]
		else:
		ch_ = ch
		ch_ = list(set(ch))
		location = "S{:d}".format(state_)
		for c in ch_:
		location=location+"_C{:d}".format(c)
		location = location+"_chan_void_vap"
		location = location+"_"+dataset
		if not figname:
		figname = location+".png"
		zs = {}
		voids = {}
		labels = []
		for c in ch_:
		(zs[c], voids[c]) = me.h5obj.getChanMeshAndData('chan_void_vap', state_, c)
		(zs[c], voids[c]) = me.h5obj.getChanMeshAndData(dataset, state_, c)
		labels.append("ch {:4d}".format(c))
		me._plotAxial(voids, zs, "Void [-]", figname, labels, location)
		me._plotAxial(voids, zs, plotYLabel, figname, labels, location)

		def plotChAxialVoid(me, state=None, ch=None, figname=None):
		""" Makes an axial plot of void for a channel

		Args:
		state (int): State number to plot (1-based). If not present, selects state with maximum.
		ch (int): Chan number to plot (1-based). If not present, selects chan with maximum.
		figname (str): Name of the figure.
		"""
		me._chanPlotMaxWrapper('chan_void_vap', "Void [-]", state, ch, figname)

		def plotChAxialQuality(me, state=None, ch=None, figname=None):
		""" Makes an axial plot of equilibrium quality for a channel

		Args:
		state (int): State number to plot (1-based). If not present, selects state with maximum.
		ch (int): Chan number to plot (1-based). If not present, selects chan with maximum.
		figname (str): Name of the figure.
		"""
		me._chanPlotMaxWrapper('chan_equilibrium_quality', "Quality [-]", state, ch, figname)

		def plotChAxial(me, state=None, ch=None, figname=None):
		""" Makes an axial plot of equilibrium quality for a channel

		Args:
		state (int): State number to plot (1-based). If not present, selects state with maximum.
		ch (int): Chan number to plot (1-based). If not present, selects chan with maximum.
		figname (str): Name of the figure.
		"""
		me._chanPlotMaxWrapper('massflux_tot', "Mass flux [kg/m**2/s]", state, ch, figname)

SubKit/process/plotChanQuality.py

0 → 100644

+18 −0

Original line number	Diff line number	Diff line
		import argparse
		from PlotTools import PlotTools

		def main():
		parser = argparse.ArgumentParser(description="Plot the axial equilibrium quality distribution in a channel. All indices are 1-based.")
		parser.add_argument('h5name', type=str, help="HDF5 file")
		parser.add_argument('--state', type=int, nargs="?", help="State to plot. Defaults to state with max quality.")
		parser.add_argument('--ch', type=int, action='append', help="""Channel to plot. Defaults to channel with max quality.
		Can plot multiple channels (e.g. --ch 1 --ch 2 --ch 3 ...)""")
		parser.add_argument('--figname', type=str, nargs="?", help="Name of figure. Defaults to name based on location and dataset type.")
		args = parser.parse_args()

		plotter = PlotTools(args.h5name)

		plotter.plotChAxialQuality(args.state, args.ch, args.figname)

		if __name__=='__main__':
		main()

SubKit/process/plotPinRadialTemp.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -6,7 +6,7 @@ def main():
		parser.add_argument('h5name', type=str, help="HDF5 file")
		parser.add_argument('--pin', type=int, nargs="?", help="Pin to plot. Defaults to pin with max centerline temperature.")
		parser.add_argument('--level', type=int, nargs="?", help="Axial level to plot. Defaults to level with max centerline temp.")
		parser.add_argument('--state', type=int, nargs="?", help="State to plot", default=1)
		parser.add_argument('--state', type=int, nargs="?", help="State to plot")
		parser.add_argument('--figname', type=str, nargs="?", help="Name of figure")
		args = parser.parse_args()

setup.py

+1 −0

Original line number	Diff line number	Diff line
		@@ -10,6 +10,7 @@ setup(name='SubKit',
		'skplot_pin_radial_temp=SubKit.process.plotPinRadialTemp:main',
		'skplot_pin_dnbr=SubKit.process.plotPinDNBR:main',
		'skplot_chan_void=SubKit.process.plotChanVoid:main',
		'skplot_chan_quality=SubKit.process.plotChanQuality:main',
		'sksummary_dnb=SubKit.process.genDNBSummary:main',]
		},
		license='MIT',