Add a lot of methods to the hdf5 interface (60dbb229) · Commits · Phathanapirom, Birdy / SubKit

MANIFEST.in

+1 −1

Original line number	Diff line number	Diff line
		include SMBuilder/input.sch
		include SMBuilder/wasp/*
		recursive-include SMBuilder/wasp *

SMBuilder/InpParse.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -6,7 +6,7 @@ import os
		import sys
		from subprocess import Popen, PIPE, check_call
		mypath = os.path.abspath(os.path.dirname(__file__))
		localWaspPath = os.path.join(mypath, "../wasp/build/install")
		localWaspPath = os.path.join(mypath, "wasp/build/install")
		CIpath = "/Users/rsk/SMBuilder/SMBuilder/wasp/build/install"
		if os.path.exists(localWaspPath):
		# This is where users should have built and installed WASP

SMBuilder/UnitConversions.py

+15 −2

Original line number	Diff line number	Diff line
		@@ -20,19 +20,32 @@ t_cm_m = 0.01 # cm/m
		t_in_m = t_in_cm*t_cm_m
		t_m_in = 1.0/t_in_m
		t_J_kJ = 0.001 # J/kJ
		t_BTU_J = 1055.06
		t_J_BTU = 1.0/t_BTU_J
		t_MBTU_J = t_BTU_J*1e6
		t_J_MBTU = 1.0/t_MBTU_J
		t_BTU_kJ = t_BTU_J*t_J_kJ
		t_kJ_BTU = 1.0/t_BTU_kJ
		t_MBTU_kJ = t_MBTU_J*t_J_kJ
		t_kJ_MBTU = 1.0/t_MBTU_kJ
		t_gal_L = 3.78541 # gal/L
		t_L_gal = 1.0/t_gal_L
		t_L_m3 = 1e-3 # L/m**3
		t_m3_L = 1.0/t_L_m3
		t_min_s = 60.0
		t_s_min = 1.0/t_min_s
		t_hr_s = 3600.0
		t_s_hr = 1.0/t_hr_s
		t_ft_m = 0.3048
		t_m_ft = 1.0/t_ft_m
		t_lbm_kg = 0.453592
		t_kg_lbm = 1.0/t_lbm_kg
		t_kg_Mlbm = t_kg_lbm/1e6
		t_Mlbm_kg = 1.0/t_kg_Mlbm
		t_bar_MPa = 0.1
		t_MPa_bar = 1.0/t_bar_MPa
		t_K_F = 1.8
		t_bar_psi = 14.5038

		def C2K(T):
		return T+273.15

SMBuilder/postProcess/CtfHDF5.py

+426 −10

Original line number	Diff line number	Diff line
		@@ -31,14 +31,200 @@ class CtfHDF5(object):
		else:
		me.numPin = me.numPin+1

		me.numCh = 1
		while True:
		if 'STATE_0001/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 """
		return me.numPin

		def getNumChan(me):
		""" Returns the total number of channels in the model """
		return me.numCh

		def isTransient(me):
		""" Returns True if this simulation is a transient """
		if 'time' in me.h5['STATE_0001']:
		return True
		else:
		return False

		def getTransientTime(me, state):
		""" Returns the transient time for this state

		Args:
		state (int) : State number (1-based)

		"""
		assert(me.isTransient())
		return me.h5['STATE_{:04d}/time'.format(state)][()][0]

		def getNumLev(me, ch):
		""" Returns the number of levels in passed channel. This includes ghost levels."""
		assert(ch>=1 and ch<=me.getNumChan())
		path = "CORE/mesh/CHAN_{:05d}/scalar_bounds".format(ch)
		return len(me.h5[path])+1

		def getPinNumLev(me, pin):
		""" Returns the number of levels in the pin

		Args:
		pin (int) : Pin index (1-based)

		"""
		return len(me.h5['CORE/mesh/PIN_{:05d}/axial_nodes'.format(pin)])

		def allChanSameLevels(me):
		""" Returns True if all channels in the model have the same number of levels."""
		nlev = None
		for ch in range(1, me.getNumChan()+1):
		path = "CORE/mesh/CHAN_{:05d}/scalar_bounds".format(ch)
		if not nlev:
		nlev = len(me.h5[path])
		else:
		thisLev = len(me.h5[path])
		if nlev!=thisLev:
		return False
		return True

		def getNumState(me):
		""" Returns the total number of states in the model """
		return me.numState

		def getNumSections(me):
		""" Returns the total number of axial sections in the model """
		return max(me.h5['/CORE/ch_sections'][()])

		def getChansInSection(me, section=1):
		""" Returns a list of the channel IDs in the passed section.

		Args:
		section (int) : Section index (1-based)

		"""
		res = []
		sections = me.h5['/CORE/ch_sections'][()]
		for i, sec in enumerate(sections):
		if section==sec:
		res.append(i+1)
		return res

		def getCrossSectionalArea(me, section=1):
		""" Returns the cross-sectional area of passed section index.

		Args:
		section (int) : Section index (1-based)

		"""
		assert(isinstance(section, int))
		area = 0.0
		chans = me.getChansInSection(section)
		return np.sum(me.h5['CORE/ch_area'][()][np.array(chans)-1])

		def _areaWeightAverage(me, dset, chans, level, state):
		""" Calculates an area-weighted average of the specified dataset

		Args:
		dset (str) : Name of the dataset
		chans (list) : List of 1-based indices of channels to include in the average
		level (int) : Axial level (0-based) at which to do the average
		state (int) : State number (1-based) to do the average

		"""
		val = 0.0
		area = 0.0
		for ch in chans:
		path = "STATE_{:04d}/CHAN_{:05d}/{:s}".format(state, ch, dset)
		areaPath = "/CORE/ch_area".format(ch)
		val = val+me.h5[path][()][level]*me.h5[areaPath][()][ch-1]
		area = area+me.h5[areaPath][ch-1]
		return val/area

		def getAvePressure(me, section=1, level=None, state=1):
		""" Returns the average pressure (area weighted)

		Args:
		section (int) : Section index from which pressure will be calculated. Defaults to 1. (1-based)
		level (int) : Axial level index in section. Defaults to outlet. (0-based)
		state (int) :: The state from which to retrieve the average pressure. (1-based)

		"""
		if not me.allChanSameLevels() or section!=1:
		raise NotImplementedError("Does not support multi-section models yet")
		if level is None:
		_level = me.getNumLev(1)-1
		else:
		_level = level
		chans = me.getChansInSection(section)
		return me._areaWeightAverage('pressure', chans, _level, state)

		def _sumChDataset(me, dset, state, chans, level):
		""" Gets the sum of the dataset for passed state, chans, and level"""
		assert(isinstance(state, int))
		assert(isinstance(chans, list))
		assert(isinstance(level, int))
		res = 0.0
		for ch in chans:
		path = "STATE_{:04d}/CHAN_{:05d}/{:s}".format(state, ch, dset)
		res = res+me.h5[path][()][level]
		return res


		def getAveMassFlux(me, section=1, level=None, state=1):
		""" Returns the average pressure (area weighted)

		Args:
		section (int) : Section index from which pressure will be calculated. Defaults to 1. (1-based)
		level (int) : Axial level index in section. Defaults to outlet. (0-based)
		state (int) :: The state from which to retrieve the average pressure. (1-based)

		"""
		if not me.allChanSameLevels() or section!=1:
		raise NotImplementedError("Does not support multi-section models yet")
		if level is None:
		_level = me.getNumLev(1)-2
		else:
		_level = level
		chans = me.getChansInSection(section)
		mdotTot = me._sumChDataset('mdot_liq', state, chans, _level)+\
		me._sumChDataset('mdot_vap', state, chans, _level)+\
		me._sumChDataset('mdot_drp', state, chans, _level)
		return mdotTot/me.getCrossSectionalArea(section)

		def getAveTemp(me, section=1, level=None, state=1):
		""" Returns the average temperature (flow weighted)

		Args:
		section (int) : Section index from which data will be calculated. Defaults to 1. (1-based)
		level (int) : Axial level index in section. Defaults to outlet. (0-based)
		state (int) :: The state from which to retrieve the data. (1-based)

		"""
		if not me.allChanSameLevels() or section!=1:
		raise NotImplementedError("Does not support multi-section models yet")
		if level is None:
		_level = me.getNumLev(1)-2
		else:
		_level = level
		temp = 0.0
		flow = 0.0
		for ch in range(1, me.getNumChan()+1):
		path = "STATE_{:04d}/CHAN_{:05d}/temperature".format(state, ch)
		flowPath = "/STATE_{:04d}/CHAN_{:05d}/".format(state, ch)
		mdotl = me.h5[flowPath+'mdot_liq'][()][_level]
		mdotv = me.h5[flowPath+'mdot_vap'][()][_level]
		mdote = me.h5[flowPath+'mdot_drp'][()][_level]
		mdotTot = mdotl+mdotv+mdote
		temp = temp +me.h5[path][()][_level]*mdotTot
		flow = flow+mdotTot
		return temp/flow


		def getPinWithMaxCenterlineTemp(me, state=None):
		""" Finds the pin with the maximum centerline temperature.

		@@ -47,9 +233,9 @@ class CtfHDF5(object):
		By default, it will search all states.

		Returns:
		state (int): State where max temp occurs
		pin (int): Pin where max temp occurs
		level (int): Level in pin where max temp occurs
		state (int): State where max temp occurs (1-based)
		pin (int): Pin where max temp occurs (1-based)
		level (int): Level in pin where max temp occurs (0-based)
		centerline (float): Value of maximum centerline temperature

		"""
		@@ -138,6 +324,44 @@ class CtfHDF5(object):

		return radialMesh, radial

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

		Args:
		pin (int) : Pin index (1-based)
		level (int) : Level index in the pin (0-based)

		"""
		if me.getNumSections()==1:
		if level==0 or level==1:
		return 1
		elif level==me.getPinNumLev(pin)-2:
		return level-1
		else:
		return level
		else:
		raise NotImplementedError("Not functioning for multi-section models currently")

		def getPinSection(me, pin, level):
		""" Returns the axial section the pin level exists in

		Args:
		pin (int) : Pin index (1-based)
		level (int) : Pin axial level (0-based)

		"""
		if me.getNumSections()==1:
		return 1
		else:
		raise NotImplementedError("Not functioning for multi-section models currently")


		def getConnChan(me, pin, level, theta):
		""" Returns the index of the channel connected to passed pin/theta index"""
		path = "CORE/mesh/PIN_{:05d}/chan_index".format(pin)
		section = me.getPinSection(pin, level)
		return me.h5[path][section-1, theta]

		def getPinWithMinDNBR(me, state=None):
		""" Returns state, pin, theta index, level index experiencing minimum DNBR in the model

		@@ -174,7 +398,17 @@ class CtfHDF5(object):
		minState = state
		return minState, minPin, minInd[1], minInd[0], minDNBR

		def getPinAxialDNBR(me, state=None, pin=None, theta=None):
		def _getPinSurface(me, dset, state, pin, theta, local):
		""" Returns axial mesh and data for passed state/pin/theta"""

		path = 'STATE_{:04d}/PIN_{:05d}/{:s}'.format(state, pin, dset)

		axial = me.h5[path+'_axial'][()][:]
		if not local:
		axial = list(np.array(axial)+me._getPinStart(pin))
		return axial, me.h5[path][()][:, theta]

		def getPinAxialDNBR(me, state=None, pin=None, theta=None, local=False):
		""" Returns the axial vector mesh and axial DNBR distribution on that mesh for selected pin

		Args:
		@@ -184,6 +418,9 @@ class CtfHDF5(object):
		was experienced.
		theta (int): The azimuthal index (0-base). If not provided, picks the minimum value in the
		azimuthal direction.
		local (bool): If True, the local pin mesh will be returned (0.0 at the rod bottom) and if
		False, the global mesh will be returned (bottom will be >0.0 if rod is not in the first
		section).

		Returns:
		axialMesh (list): Axial mesh list
		@@ -200,17 +437,33 @@ class CtfHDF5(object):
		if theta is not None:
		_theta = theta

		path = 'STATE_{:04d}/PIN_{:05d}/dnbr'.format(_state, _pin)
		return me._getPinSurface('dnbr', _state, _pin, _theta, local)

		def getPinAxialHeatFlux(me, state, pin, theta, local=False):
		""" Returns the axial mesh and axial heat flux distribution on that mesh for selected pin

		Args:
		state (int): State number to find pin (1-based).
		pin (int): Pin index to plot (1-based).
		theta (int): The azimuthal index (0-base).
		local (bool): If True, the local pin mesh will be returned (0.0 at the rod bottom) and if
		False, the global mesh will be returned (bottom will be >0.0 if rod is not in the first
		section).

		return me.h5[path+'_axial'][()][:], me.h5[path][()][:, _theta]
		Returns:
		axialMesh (list): Axial mesh list
		q (list): Heat flux distribution
		"""
		return me._getPinSurface('pin_q_tot', state, pin, theta, local)

		def _getChAxial(me, dset, state, ch):
		def _getChanAxial(me, dset, state, ch, local=False):
		""" Returns the channel dataset and axial mesh

		Args:
		dset (str) : Name of the dataset that will be retrieved
		state (int) : State number from which to retrieve data
		ch (int) : Channel ID from which to retrieve data
		local (bool): Set to True to return the local axial mesh

		Returns:
		chAxial (list): List of axial locations from which data is available
		@@ -219,19 +472,182 @@ class CtfHDF5(object):
		"""

		path = 'STATE_{:04d}/CHAN_{:05d}/{:s}'.format(state, ch, dset)
		return me.h5[path+'_mesh'][()][:], me.h5[path][()][:]
		axial = me.h5[path+'_mesh'][()][:]
		if not local:
		axial = list(np.array(axial)+me._getChanStart(ch))
		return axial, me.h5[path][()][:]

		def getChLiqEnth(me, ch, state=1):
		def getChanLiqEnth(me, ch, state=1, local=False):
		""" Returns the channel liquid enthalpy distribution and axial mesh

		Args:
		ch (int) : Channel ID from which to retrieve data
		state (int) : State number from which to retrieve data (will default to 1)
		local (bool) : If set to True, return the channel local axial mesh. Default is to return
		the global (e.g., channels in section>1 will not start at zero)

		Returns:
		chAxial (list): List of axial locations from which data is available
		chData (list): List of solution data

		"""
		return me._getChanAxial('enthalpy_liq', state, ch, local)

		def getChanEquilibriumQuality(me, ch, state=1, local=False):
		""" Returns the channel equilibrium quality distribution and axial mesh

		Args:
		ch (int) : Channel ID from which to retrieve data
		state (int) : State number from which to retrieve data (will default to 1)
		local (bool) : If set to True, return the channel local axial mesh. Default is to return
		the global (e.g., channels in section>1 will not start at zero)

		Returns:
		chAxial (list): List of axial locations from which data is available
		chData (list): List of solution data

		"""
		return me._getChanAxial('equilibrium_quality', state, ch, local)

		def getMaxChData(me, dset, state=None, ch=None):
		""" Finds max location of channel data in model

		Args:
		dset (str) : The name of the dataset (in HDF5 file) to search
		state (int) : The state index (if None, search all states)
		ch (int) : Channel index. If None search all channels in state. If state was None, search
		all channels in model.

		Returns:
		maxState (int) : The state where the value was the maximum
		maxCh (int) : The channel index where the value was the maximum
		maxLevel (int) : The level index (0-based) where the value was the maximum
		maxValue (float) : The value of the maximum

		"""

		maxVal = 0.0
		maxState = None
		maxCh = None
		maxLevel = None
		if state is not None:
		searchStates = [state]
		else:
		searchStates = list(range(1, me.getNumState()+1))

		for st in searchStates:
		for ch in range(1, me.getNumChan()+1):
		path = 'STATE_{:04d}/CHAN_{:05d}/{:s}'.format(st, ch, dset)
		assert(path in me.h5)
		maxLocal = np.max(me.h5[path][()])
		if maxLocal > maxVal:
		maxState = st
		maxVal = maxLocal
		maxCh = ch
		maxLevel = np.argmax(me.h5[path][()])

		return maxState, maxCh, maxLevel, maxVal

		def getChanVoid(me, ch, state=1, local=False):
		""" Returns the channel void distribution and axial mesh

		Args:
		ch (int) : Channel ID from which to retrieve data
		state (int) : State number from which to retrieve data (will default to 1)
		local (bool) : If set to True, return the channel local axial mesh. Default is to return
		the global (e.g., channels in section>1 will not start at zero)

		Returns:
		chAxial (list): List of axial locations from which data is available
		chData (list): List of solution data

		"""
		return me._getChanAxial('void_vap', state, ch, local)

		def getChanHeight(me, ch):
		""" Returns the axial height of the passed channel index

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

		"""
		return np.sum(me.h5['CORE/mesh/CHAN_{:05d}/scalar_heights'.format(ch)][()])

		def _getPinStart(me, pin):
		""" Gets the global axial location at the bottom of the pin

		Pins can exist in different sections and do not necessarily start in the first section.
		This procedure returns the global location at the bottom of the pin assuming section 1 starts
		at 0.0.

		Args:
		pin (int) : Pin index (1-based)

		"""
		pinConn = me.h5['CORE/mesh/PIN_{:05d}/chan_index'.format(pin)]
		zStart = 0.0
		for section in range(pinConn.shape[0]):
		if np.all(pinConn[()][section, :]==0):
		chans = me.getChansInSection(section+1)
		zStart = zStart+me.getChanHeight(chans[0])
		else:
		break
		return zStart

		def _getChanStart(me, ch):
		""" Gets the global axial location at the bottom of the channel

		Channel can exist in different sections.
		This procedure returns the global location at the bottom of the channel assuming section 1 starts
		at 0.0.

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

		"""
		mySection = me.h5['/CORE/ch_sections'][()][ch-1]
		zStart = 0.0
		for section in range(mySection-1):
		chans = me.getChansInSection(section+1)
		zStart = zStart+me.getChanHeight(chans[0])
		return zStart


		def getPressure(me, ch, state=1, local=False):
		""" Returns the channel pressure distribution and axial mesh

		Args:
		ch (int) : Channel ID from which to retrieve data.
		state (int) : State number from which to retrieve data (will default to 1)
		local (bool) : If set to True, return the channel local axial mesh. Default is to return
		the global (e.g., channels in section>1 will not start at zero)

		Returns:
		chAxial (list): List of axial locations from which data is available
		chData (list): List of solution data

		"""
		return me._getChanAxial('pressure', state, ch, local)

		def getMassFluxTot(me, ch, state=1, local=False):
		""" Returns the channel total mass flux distribution and axial mesh

		Args:
		ch (int) : Channel ID from which to retrieve data.
		state (int) : State number from which to retrieve data (will default to 1)
		local (bool) : If set to True, return the channel local axial mesh. Default is to return
		the global (e.g., channels in section>1 will not start at zero)

		Returns:
		chAxial (list): List of axial locations from which data is available
		chData (list): List of solution data

		"""
		return me._getChAxial('enthalpy_liq', state, ch)
		axial, mdotLiq = me._getChanAxial('mdot_liq', state, ch)
		axial, mdotVap = me._getChanAxial('mdot_vap', state, ch)
		axial, mdotDrp = me._getChanAxial('mdot_drp', state, ch)
		if not local:
		axial = list(np.array(axial)+me._getChanStart(ch))
		mdotTot = np.array(mdotLiq)+np.array(mdotVap)+np.array(mdotDrp)
		return axial, list(mdotTot/me.h5['CORE/ch_area'][()][ch-1])

SMBuilder/postProcess/CtfPlotTools.py

+31 −0

Original line number	Diff line number	Diff line
		@@ -120,3 +120,34 @@ class CtfPlotTools(object):
		#ax.legend()
		plt.savefig(figname_)
		plt.close(fig)

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

		def _plotAxial(me, x, y, xname, figname, title=None, yname="Axial [m]"):
		""" Makes an axial plot"""
		fig, ax = plt.subplots()
		if title:
		plt.title(title)
		ax.grid()
		ax.set_ylabel(yname)
		ax.set_xlabel(xname)
		ax.plot(x, y)
		plt.savefig(figname)
		plt.close(fig)

		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.
		"""
		stateMax, chMax, levMax, voidMax = me.h5obj.getMaxChData('void_vap', state, ch)
		(z, void) = me.h5obj.getChanVoid(chMax, stateMax)
		location = me._chName(stateMax, chMax)
		if not figname:
		figname = location+"_void_vap.png"
		me._plotAxial(void, z, "Void [-]", figname, location)

Admin message