Add script for generating time-dependent channel data (a7316f05) · Commits · Phathanapirom, Birdy / SubKit

SubKit/process/Hdf5Tools.py

+89 −12

Original line number	Diff line number	Diff line
		@@ -42,7 +42,7 @@ class Hdf5Interface(object):
		me._transientTimes = []
		if me._isTransient:
		for state in range(me.numState):
		me._transientTimes.append(me.h5['STATE_{:04d}/time'.format(state+1)][()])
		me._transientTimes.append(me.h5['STATE_{:04d}/time'.format(state+1)][0])

		# Takes state number (1-based) and returns a list of dataset names
		me._stateDatasets = {}
		@@ -131,6 +131,7 @@ class Hdf5Interface(object):
		if 'mesh_type' in dataset.attrs:
		me.chanDatasetMeshTypes[datasetName] = dataset.attrs['mesh_type'][0]


		def getNumPin(me):
		""" Returns the total number of pins in the model """
		return me.numPin
		@@ -496,7 +497,7 @@ class Hdf5Interface(object):
		else:
		return "scalar"

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

		Returns a vector of axial data.
		@@ -505,12 +506,17 @@ class Hdf5Interface(object):
		dataset (str) : Name of the dataset from which data will be extracted
		state (int) : 1-based State index from which data will be extracted
		chan (int) : 1-based channel index
		level (int) : 0-based level index. If not provided, return a vector with data for all levels.
		"""
		me._stateValid(state)
		me._datasetValid(dataset)
		me._chanIndexValid(chan)
		numLevel = me.getChanDatasetNumLev(chan, dataset)
		if level is None:
		return me.h5['/STATE_{:04d}/{:s}'.format(state, dataset)][list(range(numLevel)), chan-1][()]
		else:
		me._chLevelValid(dataset, state, chan, level)
		return me.h5['/STATE_{:04d}/{:s}'.format(state, dataset)][()][level, chan-1]

		def getPinAxialDataset(me, dataset, state, pin):
		""" Retrieve pin dataset from the h5 file
		@@ -1180,12 +1186,64 @@ class Hdf5Interface(object):
		me._chanIndexValid(chan)
		return me.getChanMeshAndData(dataset, state, chan)

		def getMaxChanData(me, dataset, state=None):
		def datasetInState(me, dataset, state):
		""" Returns True if the dataset exists in the passed state index

		Args:
		dataset (str) : Name of the dataset
		state (int) : 1-based state index

		"""
		me._stateValid(state)
		return dataset in me._stateDatasets[state]


		def getDatasetNames(me, state):
		""" Returns the dataset names in the passed state index

		Args:
		state (int) : 1-based state index
		"""
		me._stateValid(state)
		return me._stateDatasets[state]

		def datasetInAllStates(me, dataset):
		""" Returns True if the dataset name exists in all states in the simulation

		Returns False if it exists in some, but not all.
		Returns False if it exists in none.

		Args:
		dataset (str) : A valid dataset name in the hdf5 file

		"""
		for state in me._stateDatasets:
		if dataset not in me._stateDatasets[state]:
		return False
		return True

		def datasetInAnyState(me, dataset):
		""" Returns True if the dataset name exists in any state in the simulation

		Returns False if it exists in none

		Args:
		dataset (str) : Dataset name to search for in the hdf5 file

		"""
		for state in me._stateDatasets:
		if dataset in me._stateDatasets[state]:
		return True
		return False

		def getMaxChanData(me, dataset, state=None, ch=None, level=None):
		""" Finds max location of channel data in model

		Args:
		dataset (str) : The name of the dataset (in HDF5 file) to search
		state (int) : The state index (if None, search all states)
		ch (int) : The 1-based channel index to search for the max (None means search all)
		level (int) : The 0-based level index to search for the max (None means search all)

		Returns:
		maxState (int) : The state where the value was the maximum
		@@ -1204,17 +1262,36 @@ class Hdf5Interface(object):
		else:
		searchStates = list(range(1, me.getNumState()+1))

		if ch is None:
		searchChans = list(range(me.getNumChan()))
		else:
		assert(isinstance(ch, int))
		searchChans = ch-1

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

		for st in searchStates:
		ma = np.ma.array(me.h5['STATE_{:04d}/{:s}'.format(st, dataset)], mask=mask)
		if level is None:
		ma = np.ma.array(me.h5['STATE_{:04d}/{:s}'.format(st, dataset)][()][:, searchChans], mask=mask[:, searchChans])
		else:
		ma = np.ma.array(me.h5['STATE_{:04d}/{:s}'.format(st, dataset)][()][level, searchChans], mask=mask[level, searchChans])
		maxval = ma.max()
		if maxVal is None or maxval>maxVal:
		maxVal = maxval
		maxState = st
		if ch is not None and level is not None:
		maxChan = ch-1
		maxLevel = level
		elif ch is not None:
		maxLevel = ma.argmax()
		maxChan = ch-1
		elif level is not None:
		maxChan = ma.argmax()
		maxLevel = level
		else:
		(maxLevel, maxChan) = np.unravel_index(ma.argmax(), ma.shape)

		return maxState, maxChan+1, maxLevel, maxVal

SubKit/process/SummaryTools.py

+96 −7

Original line number	Diff line number	Diff line
		@@ -3,8 +3,9 @@ import sys
		import os
		sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '../../')))
		import SubKit
		import SubKit.utils.UnitConversions as unit
		from SubKit.utils.UnitConversions import *
		from SubKit.utils.Timer import eTime
		import warnings

		class SummaryTools(object):
		""" Set of tools for generating text summary files from the HDF5 file
		@@ -18,6 +19,94 @@ class SummaryTools(object):
		print me.time.elapsed(), "Initializing file"
		me.h5 = Hdf5Interface(h5file)

		# Units expected in HDF5 file and their conversion constants
		me.unitConversions = {'kJ/kg' : [t_kJ_BTU/t_kg_lbm, 'BTU/lbm' ],
		'n/a' : [1, '---' ],
		'kg/s' : [t_kg_lbm/t_s_hr, 'lbm/hr' ],
		'bar' : [t_bar_psi, 'psi' ],
		'kg/m3' : [t_kg_lbm/t_m_ft3, 'lbm/ft**3' ],
		'C' : [None, 'F' ],
		'm/s' : [t_m_ft, 'ft/s' ],
		'kW/m2/K' : [t_kJ_MBTU/t_m_ft2/t_K_F/t_s_hr, 'MBTU/hr/ft**2/F'],
		'm' : [t_m_ft, 'ft' ],
		'kW/m2' : [t_kJ_MBTU/t_m_ft2/t_s_hr, 'MBTU/hr/ft**2' ],
		'J/kg' : [t_J_BTU/t_kg_lbm, 'BTU/lbm' ]}

		def genChanTimeSummary(me, dataset, fname=None, ch=None, level=None, units='si'):
		""" Prints summary information for a channel dataset

		Args:
		dataset (str) : Name of the dataset (corresponds to name in hdf5 file).
		fname (str) : The name of the file (defaults to <dataset>_CX_LX_summary.txt) if not provided where C and L are
		channel (1-based) and level (0-based) indices
		ch (int) : The (1-based) channel index to print. If not provided, finds max value in the dataset in the model
		(space and time).
		level (int) : The (1-based) level index to print. If not provided, finds the max value in the dataset in
		the model (space and time).
		units (str) : Can be 'si' or 'british'

		"""
		if not me.h5.datasetInAnyState(dataset):
		raise RuntimeError("Dataset: "+dataset+" not found in any STATE in HDF5 file.")
		if not dataset.startswith('chan_'):
		raise RuntimeError("Dataset: "+dataset+" is not a valid channel dataset.")
		if level is None:
		zeroLevel = level
		else:
		zeroLevel = level-1
		_stateMax, _chMax, _levelMax, _valMax = me.h5.getMaxChanData(dataset, ch=ch, level=zeroLevel)
		if not fname:
		_fname = '{:s}_C{:d}_L{:d}_summary.txt'.format(dataset, _chMax, _levelMax+1)
		else:
		_fname = fname
		f = open(_fname, 'w')
		f.write('#========================================================\n')
		f.write('# Channel with max: {:8d}\n'.format(_chMax))
		f.write('# Level with max: {:8d}\n'.format(_levelMax+1))
		f.write('# Dataset: {:s}\n'.format(dataset))
		f.write('#--------------------------------------------------------\n')
		if me.h5.isTransient():
		time = 'Time [s]'
		else:
		time = 'State '
		_units = me.h5.getDatasetUnits(dataset)
		if units!='si' and units!='british':
		raise RuntimeError("Unexpected units: "+units+", must be 'si' or 'british'")
		doConversion = False
		if units=='british':
		if _units not in me.unitConversions:
		warnings.warn('Unit type of dataset: '+dataset+' is not an expected type: '+_units+' so not converting to british')
		else:
		conversion = me.unitConversions[_units][0]
		if _units=='C':
		dataIsTemp = True
		else:
		dataIsTemp = False
		_units = me.unitConversions[_units][1]
		doConversion = True

		f.write('#{:>15s}{:>25s}\n'.format(time, dataset+' ['+_units+']'))
		for state in range(1, me.h5.getNumState()+1):
		if me.h5.datasetInState(dataset, state):
		val = me.h5.getChanDataset(dataset, state, _chMax, _levelMax)
		if units=='british' and doConversion:
		if dataIsTemp:
		val = C2F(val)
		else:
		val = val*conversion
		if me.h5.isTransient():
		if isinstance(val, int):
		f.write('{:15.5e}{:25d}\n'.format(me.h5.getTransientTime(state), val))
		else:
		f.write('{:15.5e}{:25.5e}\n'.format(me.h5.getTransientTime(state), val))
		else:
		if isinstance(val, int):
		f.write('{:15d}{:25d}\n'.format(state, val))
		else:
		f.write('{:15d}{:25.5e}\n'.format(state, val))
		f.close()


		def genDNBSummary(me, fname=None, units=None):
		""" Prints summary information for limiting rod/channel"""

		@@ -58,28 +147,28 @@ class SummaryTools(object):
		axial, dnbr = me.h5.getPinAxialDNBR(state, minPin, theta)
		sysPressure = me.h5.getAvePressure(state=state)
		if convert:
		sysPressure = sysPressure*unit.t_bar_psi
		sysPressure = sysPressure*t_bar_psi
		inletTemp = me.h5.getAveTemp(state=state)
		if convert:
		inletTemp = unit.C2F(inletTemp)
		inletTemp = C2F(inletTemp)
		inletMflux = me.h5.getAveMassFlux(state=state)
		if convert:
		inletMflux = inletMfluxunit.t_kg_lbm/unit.t_m_ft*2/unit.t_s_hr
		inletMflux = inletMfluxt_kg_lbm/t_m_ft*2/t_s_hr
		if convert:
		location = axial[level]*unit.t_m_in
		location = axial[level]*t_m_in
		else:
		location = axial[level]
		axial, mflux = me.h5.getMassFluxTot(hotChID, state)
		chLevel = me.h5.getChanLevelFromPinLevel(minPin, level, hotChID)
		mfluxHot = mflux[chLevel]
		if convert:
		mfluxHot = mfluxHotunit.t_kg_Mlbm/unit.t_s_hr/unit.t_m_ft*2
		mfluxHot = mfluxHott_kg_Mlbm/t_s_hr/t_m_ft*2
		axial, quality = me.h5.getChanEquilibriumQuality(hotChID, state)
		xHot = quality[chLevel]
		axial, q = me.h5.getPinAxialHeatFlux(state, minPin, theta)
		qHot = max(q[level], 0.0)
		if convert:
		qHot = qHotunit.t_kJ_MBTU/unit.t_s_hr/unit.t_m_ft*2
		qHot = qHott_kJ_MBTU/t_s_hr/t_m_ft*2
		if me.h5.isTransient():
		time = me.h5.getTransientTime(state)
		else:

SubKit/process/genChanTimeDataSummary.py

0 → 100644

+19 −0

Original line number	Diff line number	Diff line
		import argparse
		from SummaryTools import SummaryTools

		def main():
		parser = argparse.ArgumentParser(description="Generates a summary of time-dependent channel solution data.")
		parser.add_argument('h5name', type=str, help="HDF5 file")
		parser.add_argument('dataset', type=str, help="Valid dataset name in HDF5 file")
		parser.add_argument('--ch', type=int, help="Channel index (1-based) to plot (if not provided defaults to channel with max dataset value.")
		parser.add_argument('--level', type=int, help="Level index (0-based) to plot (if not provided defaults to level with max dataset value.")
		parser.add_argument('--units', choices=["si", "british"], type=str, help="Selector for physical units (default si)")
		parser.add_argument('--fname', type=str, help="Output summary file name")
		args = parser.parse_args()

		s = SummaryTools(args.h5name)

		s.genChanTimeSummary(args.dataset, fname=args.fname, ch=args.ch, level=args.level, units=args.units)

		if __name__=="__main__":
		main()

setup.py

+2 −1

Original line number	Diff line number	Diff line
		@@ -17,7 +17,8 @@ setup(name='SubKit',
		'skplot_pin_temp_time=SubKit.process.plotPinTempTime:main',
		'skplot_pin_surf_temp_axial=SubKit.process.plotPinTempAxial:main',
		'sk_gen_from_template=SubKit.utils.gen_from_template:main',
		'sksummary_dnb=SubKit.process.genDNBSummary:main',]
		'sksummary_dnb=SubKit.process.genDNBSummary:main',
		'sksummary_chan_time_data=SubKit.process.genChanTimeSummaryData:main']
		},
		license='MIT',
		packages=find_packages(),

tests/unitTests/test_Hdf5Tools.py

+80 −0

Original line number	Diff line number	Diff line
		@@ -104,6 +104,51 @@ class test_Hdf5Tools(unittest.TestCase):
		for i in range(len(sec2Mom)):
		me.assertAlmostEqual(z[i], sec2Mom[i])

		def test_getChanDataset(me):

		# Get vector data
		temp = me.obj.getChanDataset('chan_temp', 2, 1)
		gold = [292.7112419154088, 294.1329178774784, 295.2869793451498, 296.2703916610701, 297.1424651121381]
		for i, t in enumerate(temp):
		me.assertAlmostEqual(t, gold[i])

		# Get scalar data
		level = 2
		me.assertAlmostEqual(me.obj.getChanDataset('chan_temp', 2, 1, level), gold[level])

		def test_dataset_checkers(me):
		# Tests datasetInState, getDatasetNames

		# Datasets that can be found in files
		datasets = ['chan_enthalpy_liq', 'chan_enthalpy_vap', 'chan_equilibrium_quality', 'chan_flow_regime',
		'chan_liq_subcooling', 'chan_mdot_drp', 'chan_mdot_liq', 'chan_mdot_vap', 'chan_pressure',
		'chan_rho_liq', 'chan_rho_vap', 'chan_temp', 'chan_velocity_drp', 'chan_velocity_liq',
		'chan_velocity_vap', 'chan_void_vap', 'pin_axial_centerline_temp', 'pin_axial_gap_conductance',
		'pin_axial_gap_contact_pressure', 'pin_axial_gap_thickness', 'pin_surf_dnbr', 'pin_surf_heat_transfer_mode',
		'pin_surf_htc_liq', 'pin_surf_q_tot', 'pin_surf_steaming', 'pin_surf_tke', 'pin_temp']

		numState = 2

		for state in range(1, numState+1):
		fileDatasets = me.obj.getDatasetNames(state)
		me.assertEqual(len(fileDatasets), len(datasets))
		for test in fileDatasets:
		me.assertTrue(test in datasets)


		# Currently datasets exist in all states or none, so cannot check if in one state and not another
		# We may add this feature in the future.
		for dataset in datasets:
		me.assertTrue(me.obj.datasetInAllStates(dataset))
		me.assertTrue(me.obj.datasetInAnyState(dataset))
		for state in range(1, numState+1):
		me.assertTrue(me.obj.datasetInState(dataset, state))


		me.assertFalse(me.obj.datasetInAllStates('junkname'))
		me.assertFalse(me.obj.datasetInAnyState('junkname'))
		me.assertFalse(me.obj.datasetInState('junkname', 1))

		def test_getChanDatasetNumLev(me):

		# Datasets on scalar mesh
		@@ -167,6 +212,41 @@ class test_Hdf5Tools(unittest.TestCase):
		me.assertEqual(level, 3)
		me.assertAlmostEqual(maxval, 0.2)

		# Only look in channel 1
		state, ch, level, maxval = me.obj.getMaxChanData('chan_temp', state=2, ch=1)
		me.assertEqual(state, 2)
		me.assertEqual(ch, 1)
		me.assertEqual(level, 4)
		me.assertAlmostEqual(maxval, 297.1424651121381)

		# Leave out the state, answer should be the same
		state, ch, level, maxval = me.obj.getMaxChanData('chan_temp', ch=1)
		me.assertEqual(state, 2)
		me.assertEqual(ch, 1)
		me.assertEqual(level, 4)
		me.assertAlmostEqual(maxval, 297.1424651121381)

		# Only look at level 1
		state, ch, level, maxval = me.obj.getMaxChanData('chan_temp', state=2, level=1)
		me.assertEqual(state, 2)
		me.assertEqual(ch, 10)
		me.assertEqual(level, 1)
		me.assertAlmostEqual(maxval, 301.9729903325481)

		# Leave out the state, answer should be the same
		state, ch, level, maxval = me.obj.getMaxChanData('chan_temp', level=1)
		me.assertEqual(state, 2)
		me.assertEqual(ch, 10)
		me.assertEqual(level, 1)
		me.assertAlmostEqual(maxval, 301.9729903325481)

		# Specify ch/level
		state, ch, level, maxval = me.obj.getMaxChanData('chan_temp', state=2, ch=2, level=1)
		me.assertEqual(state, 2)
		me.assertEqual(ch, 2)
		me.assertEqual(level, 1)
		me.assertAlmostEqual(maxval, 293.8922181191907)


		def test_getRadialTemp(me):
		# Get a radial temperature distribution without centerline