Add plot tool for pin temperature data over time (2d60a1f7) · Commits · Phathanapirom, Birdy / SubKit

SubKit/process/Hdf5Tools.py

+147 −4

Original line number	Diff line number	Diff line
		@@ -67,11 +67,12 @@ class Hdf5Interface(object):
		# 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)
		me.pinInternalMask[:, :, :, :] = True
		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
		me.pinInternalMask[0:nsec, 0:nlev, 0:nrad, pin-1] = False

		# Create a mask for pin axial data
		if 'STATE_0001/pin_axial_centerline_temp' in me.h5:
		@@ -120,6 +121,18 @@ class Hdf5Interface(object):
		assert(me.isTransient())
		return me.h5['STATE_{:04d}/time'.format(state)][()]

		def getTransientTimes(me):
		""" Returns a list of times in the transient

		Will be None if not a transient.
		"""
		if me.isTransient():
		times = []
		for state in range(1, me.getNumState()+1):
		times.append(me.getTransientTime(state))
		return times


		def _sectionValid(me, section):
		if section<1 or section>me.getNumSections():
		raise ValueError("Section number "+str(section)+" is not found in the model")
		@@ -135,6 +148,11 @@ class Hdf5Interface(object):
		if pin < 1 or pin > me.getNumPin():
		raise ValueError("Pin index "+str(pin)+" is not found in the model")

		def _pinRadialIndexValid(me, pin, radial):
		""" Used for checking validity of radial index """
		if radial<0 or radial>=me.getPinNumRadial(pin):
		raise ValueError("Pin radial index "+str(radial)+" for pin "+str(pin)+" is not in the pin")

		def _pinHasCouplingData(me, pin):
		""" Checks if the pin has coupling mesh data """
		if 'reconstruction' in me.h5['CORE/mesh/PIN_{:05d}'.format(pin)]:
		@@ -644,6 +662,42 @@ class Hdf5Interface(object):
		level (int) : The axial level (0-based) in the pin where the max occurs
		val (float) : The max value
		"""
		if state is None:
		states = list(range(1, me.getNumState()+1))
		else:
		states = [state]
		maxVal = None
		for state_ in states:
		ma = np.ma.array(me.h5['STATE_{:04d}/{:s}'.format(state_, dataset)], mask=me.pinAxialMask)
		if pin is not None:
		# A pin index was passed, so mask all other pins except this one
		for p in range(me.getNumPin()):
		if p+1!=pin:
		ma[:, p] = np.ma.masked
		maxval = ma.max()
		if maxVal is None or maxval>maxVal:
		maxVal = maxval
		maxState = state_
		(maxLevel, maxPin) = np.unravel_index(ma.argmax(), ma.shape)
		return maxState, maxPin+1, maxLevel, maxVal

		def _getPinInternalMax(me, state=None, pin=None):
		""" Use to find the max temperature in the pin_temp dataset

		Args:
		state (int) : The state number (1-based) in which to look for the max, or look in all states
		if None
		pin (int) : The pin number (1-based) in which to look for the mask, or look in all pins if None

		Returns:
		state (int) : The state number (1-based) where the max occurs (will be same as passed
		state if state was passed)
		pin (int) : The pin index (1-based) where the max occurs
		sector (int) : The sector index (0-based) where the max occurs
		rad (int) : The radial index (0-based) where the max occurs
		level (int) : The axial level (0-based) in the pin where the max occurs
		val (float) : The max value
		"""
		if state is None:
		states = list(range(1, me.getNumState()+1))
		else:
		@@ -652,15 +706,21 @@ class Hdf5Interface(object):
		pins = list(range(1, me.getNumPin()+1))
		else:
		pins = [pin]

		maxVal = None
		for state_ in states:
		ma = np.ma.array(me.h5['STATE_{:04d}/{:s}'.format(state_, dataset)], mask=me.pinAxialMask)[:, np.array(pins)-1]
		ma = np.ma.array(me.h5['STATE_{:04d}/pin_temp'.format(state_)], mask=me.pinInternalMask)
		if pin is not None:
		# A pin index was passed, so mask all other pins except this one
		for p in range(me.getNumPin()):
		if p+1!=pin:
		ma[:, :, :, p] = np.ma.masked
		maxval = ma.max()
		if maxVal is None or maxval>maxVal:
		maxVal = maxval
		maxState = state_
		(maxLevel, maxPin) = np.unravel_index(ma.argmax(), ma.shape)
		return maxState, maxPin+1, maxLevel, maxVal
		(maxSect, maxLevel, maxRad, maxPin) = np.unravel_index(ma.argmax(), ma.shape)
		return maxState, maxPin+1, maxSect, maxRad, 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
		@@ -777,6 +837,32 @@ class Hdf5Interface(object):

		return radialMesh, Tave

		def _getPinMaxTempTime(me, state=None, pin=None):
		""" Returns the location where max pin temperature occurs through time.

		Args:
		state (int) : 1-based state index. If not passed, searches all states.
		pin (int) : 1-based pin index. If not passed, searches all pins.

		Returns:
		state (int) : 1-based state where max temp occurs.
		pin (int) : 1-based pin index where max temp occurs.
		sect (int) : 0-based sector index where max temp occurs. None will be returned if the max is at the centerline.
		radial (int) : Radial index where max temp occurs. 0 means centerline.
		level (int) : 0-based level index where max temp occurs.
		temp (float) : The value of the maximum
		"""
		# Find location with max centerline
		stateCent, pinCent, levelCent, maxTcent = me.getPinWithMaxCenterlineTemp(pin=pin)
		# Find location with non-centerline max
		state, pin, sect, rad, level, maxT = me._getPinInternalMax(pin=pin)

		if maxTcent>maxT:
		return stateCent, pinCent, None, 0, levelCent, maxTcent
		else:
		return state, pin, sect, rad+1, level, maxT


		def getChanLevelFromPinLevel(me, pin, level, chan):
		""" Takes the pin index and level and returns the scalar level index in the adjacent channel

		@@ -1096,3 +1182,60 @@ class Hdf5Interface(object):
		mdotTot = np.array(mdotLiq)+np.array(mdotVap)+np.array(mdotDrp)
		return axial, list(mdotTot/me.h5['CORE/chan_area'][()][chan-1])

		def getPinAxialTimeData(me, dataset, pin, level):
		""" Returns time series data for single location in pin axial data

		Pin axial data is data that is dependent on axial location only.

		Args:
		dataset (str) : The name of the dataset in the HDF5 file.
		pin (int) : The 1-based pin index
		level (int) : The 0-based level index.

		Returns:
		time (float or int) : A list of transient times (if transient) or state indices (if depletion)
		data (float) : A list of data corresponding to the time data

		"""
		me._pinIndexValid(pin)
		me._pinLevelValid(pin, level)
		me._datasetValid(dataset)
		times = []; data = []
		for state in range(1, me.getNumState()+1):
		if me.isTransient():
		times.append(me.h5['STATE_{:04d}/time'.format(state)][()])
		else:
		times.append(state)
		data.append(me.h5['STATE_{:04d}/{:s}'.format(state, dataset)][()][level, pin-1])
		return times, data

		def getPinInternalTimeData(me, pin, theta, radial, level):
		""" Returns time series data for single location in pin internal data

		The pin internal data currently only covers pin_temp. Takes a single spatial location and
		returns a list of times and a list of values for those times.

		Args:
		pin (int) : 1-based pin index
		theta (int) : 0-based sector index
		radial (int) : 0-based radial node index
		level (int) : 0-based level index

		Returns:
		time [real/int] : If a transient, this will be a list of times in seconds and if it is not then
		it will be a list of integers (one per state)
		data [real] : List of values at times

		"""
		me._pinIndexValid(pin)
		me._pinSectorValid(pin, theta)
		me._pinRadialIndexValid(pin, radial)
		me._pinLevelValid(pin, level)
		times = []; temps = []
		for state in range(1, me.getNumState()+1):
		if me.isTransient():
		times.append(me.h5['STATE_{:04d}/time'.format(state)][()])
		else:
		times.append(state)
		temps.append(me.h5['STATE_{:04d}/pin_temp'.format(state)][()][theta, level, radial, pin-1])
		return times, temps

SubKit/process/PlotTools.py

+88 −0

Original line number	Diff line number	Diff line
		@@ -44,6 +44,7 @@ class PlotTools(object):

		return state_, pin_, level_


		def plotPinRadialTemp(me, state=None, pin=None, level=None, figname=None):
		""" Make a radial temperature plot in a pin

		@@ -78,6 +79,93 @@ class PlotTools(object):
		plt.savefig(figname_)
		plt.close(fig)

		def plotPinTimePlot(me, pin=None, theta=None, rad=None, level=None, figname=None):
		""" Makes a time-plot of pin internal data.

		Can provide a list of data (in which case all location arguments must be the same length), a single location,
		or no location at all (all must be None) to get the max temperature in space and time for plotting.

		Args:
		pin (int) : 1-based pin index. Can be scalar or list. Picks max if not provided.
		theta (int) : 0-based sector index. Can be scalar or list. Picks max if not provided.
		rad (int) : 1-based radial index. Can be scalar or list. Picks max if not provided. 0 means centerline and
		negative values count back from the surface with -1 being the surface value.
		level (int) : 1-based level index. Can be scalar or list. Picks max if not provided.
		figname (str) : Name of the figure. Picks name based on location if not provided. Note that R0 means the location
		was at the centerline and R1 means the first radial node.

		"""
		if pin is None or theta is None or rad is None or level is None:
		if pin is not None or theta is not None or rad is not None or level is not None:
		raise RuntimeError("If pin/theta/rad/level is None, then all must be None")
		# All are None, so find the max
		state, pin, theta, rad, level, temp = me.h5obj._getPinMaxTempTime()
		states = [state]; pins = [pin]; thetas = [theta]; rads = [rad]; levels = [level]
		else:
		if isinstance(pin, int):
		pins = [pin]
		else:
		pins = pin
		if isinstance(theta, int):
		thetas = [theta]
		else:
		thetas = theta
		if isinstance(rad, int):
		rads = [rad]
		else:
		rads = rad
		if isinstance(level, int):
		levels = [level]
		else:
		levels = level
		# Ensure all lists are the same length
		if len(pins)!=len(rads) or len(pins)!=len(levels) or len(pins)!=len(thetas):
		raise RuntimeError("Length of pin/theta/rad/level locations must be of equal size.")
		# Convert negative radial locations to positive
		for i, p in enumerate(pins):
		if rads[i]<0:
		nrad = me.h5obj.getPinNumRadial(p)
		rads[i] = nrad+rads[i]+1

		# Make all indices 1-based for I/O
		location = "P{:d}".format(pins[0])
		if len(pins)>1:
		for p in pins[1:]:
		location = location+"_P{:d}".format(p)
		for l in levels:
		location = location+"_L{:d}".format(l+1)
		for r in rads:
		location = location+"_R{:d}".format(r)
		for t in thetas:
		location = location+"_T{:d}".format(t+1)
		location = location+"_pin_temp"
		if not figname:
		figname = location+".png"


		fig, ax = plt.subplots()
		plt.title(location)
		ax.grid()
		if me.h5obj.isTransient():
		ax.set_xlabel("Time [s]")
		else:
		ax.set_xlabel("State [-]")
		ax.set_ylabel("Temperature [C]")
		for i, p in enumerate(pins):
		if rads[i]==0:
		times, temps = me.h5obj.getPinAxialTimeData('pin_axial_centerline_temp', p, levels[i])
		else:
		times, temps = me.h5obj.getPinInternalTimeData(p, thetas[i], rads[i]-1, levels[i])
		label = 'P{:d}_L{:d}_R{:d}_T{:d}'.format(p, levels[i]+1, rads[i], thetas[i]+1)
		ax.plot(times, temps, label=label)
		if not me.h5obj.isTransient():
		xint = range(times[0], times[-1]+1)
		ax.set_xticks(xint)
		ax.legend()
		plt.savefig(figname)
		plt.close(fig)


		def plotPinAxialDNBR(me, state=None, pin=None, theta=None, figname=None):
		""" Makes an axial plot of DNBR for a pin

SubKit/process/plotPinTempTime.py

0 → 100644

+42 −0

Original line number	Diff line number	Diff line
		import argparse
		from PlotTools import PlotTools
		import numpy as np

		def main():
		parser = argparse.ArgumentParser(description="Plot the temperature of a rod location through time. "+
		"A location will be printed for all STATE points in the results file. If this is a transient, "+
		"the x-axis will be time (s) and if it is a depletion, the x-axis will be STATE index. Multiple locations "+
		"can be plotted in the same plot. Provide no locations to plot the location with max temperature through "+
		"time. If locations are provided, all must be provided and the number of pin/rad/level/theta must be equal in size.")
		parser.add_argument('h5name', type=str, help="HDF5 file")
		parser.add_argument('--pin', type=int, action='append', help="Pin to plot. Can provide multiple to plot multiple pins.")
		parser.add_argument('--rad', type=int, action='append', help="Radial location to plot. Can provide multiple to plot "+
		"multiple radial locations. Enter '0' to plot the centerline temperature. Can also specify negative number to "+
		"count backwards from the surface (e.g., -1 would plot the outer clad surface temperature for a fuel rod and "+
		"-2 would be the inner clad surface temperature).")
		parser.add_argument('--theta', type=int, action='append', help='The azimuthal sector index to plot. Can provide multiple.')
		parser.add_argument('--level', type=int, action='append', help="Axial level to plot. Can specify multiple axial levels.")
		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)

		if args.pin is None or args.rad is None or args.theta is None or args.level is None:
		if args.pin is not None or args.rad is not None or args.theta is not None or args.level is not None:
		raise RuntimeError("One of --pin, --rad, --theta, or --level was provided. Must pass each option if any are specified.")
		thetas = None
		levels = None
		if args.pin and args.rad and args.theta and args.level:
		if len(args.pin)!=len(args.rad):
		raise RuntimeError("Length of pin and rad lists provided do not match: {:d}/={:d}".format(len(args.pin), len(args.rad)))
		if len(args.pin)!=len(args.level):
		raise RuntimeError("Length of pin and level lists provided do not match: {:d}/={:d}".format(len(args.pin), len(args.level)))
		if len(args.pin)!=len(args.theta):
		raise RuntimeError("Length of pin and theta lists provided do not match: {:d}/={:d}".format(len(args.pin), len(args.theta)))
		thetas = list(np.array(args.theta)-1)
		levels = list(np.array(args.level)-1)

		plotter.plotPinTimePlot(args.pin, thetas, args.rad, levels, args.figname)

		if __name__=='__main__':
		main()

setup.py

+1 −0

Original line number	Diff line number	Diff line
		@@ -13,6 +13,7 @@ setup(name='SubKit',
		'skplot_chan_quality=SubKit.process.plotChanQuality:main',
		'skplot_pin_coupling_tke=SubKit.process.plotCouplingSurfTKE:main',
		'skplot_pin_coupling_temp=SubKit.process.plotCouplingSurfTemp:main',
		'skplot_pin_temp_time=SubKit.process.plotPinTempTime:main',
		'sk_gen_from_template=SubKit.utils.gen_from_template:main',
		'sksummary_dnb=SubKit.process.genDNBSummary:main',]
		},

tests/unitTests/dummy_transient.h5

0 → 100644

+127 KiB

File added.

No diff preview for this file type.

View file