Adds a plot tool for channel axial liquid temperature distribution and one for... (81114925) · Commits · Phathanapirom, Birdy / SubKit

SubKit/build/Model.py

+30 −1

Original line number	Diff line number	Diff line
		@@ -94,6 +94,11 @@ class Model:
		me.editOptions["convergence"] = True
		me.editOptions["gap_edits"] = False

		# Model maps not supplied by default
		me.rodMap = None
		me.chanMap = None
		me.symOpt = None

		# Steady state convergence criteria
		me.steadyConvergenceCriteria = {'massBalance': 0.1,
		'energyBalance': 0.1,
		@@ -150,9 +155,10 @@ class Model:
		- "if97_direct" (water) -- Use IF97 thermophysical properties for water (slower than tables).
		- "asme_1968" (water) -- Use lookup tables generated from ASME 1968 thermophysical properties for water.
		- "flibe" (molten salt) -- Use properties for FLiBe
		- "custom" -- User will provide a set of custom fluid properties in thermophysical_properties.dat file.

		"""
		validEntries = ['if97_lookup', 'if97_direct', 'asme_1968', 'flibe']
		validEntries = ['if97_lookup', 'if97_direct', 'asme_1968', 'flibe', 'custom']
		assert(fluidprops in validEntries)
		me.fluidProperties = fluidprops

		@@ -951,6 +957,29 @@ class Model:
		return secID
		assert(False) # Fail if not found in any sections

		def setCoreMaps(me, rodMap, chanMap, assemMap, symOpt=1):
		""" Sets the rod/channel layouts for the model

		If the model has a rod-bundle type of layout, a set of maps can be supplied.
		The rodMap will be a 2D array (list of lists) that takes pin row/col and returns the assembly
		index which owns the pin. The chanMap will similarly be a 2D array (list of lists) that takes chan row/col
		and returns the assembly index which owns the channel. The chanMap should be one row/col larger
		than the rod map and it is assumed that the channels will sit between the rods.
		The assemMap takes the assembly row/col and returns the assembly index.
		This information must be specified if the user wishes to write a rod VTK file, a VERA HDF5 file,
		or for CTF to be used in a coupled simulation.

		Args:
		rodMap (list) : 2D list specifying rod layout
		chanMap (list) : 2D list specifying channel layout
		assemMap (list) : 2D list specifying assembly layout
		symOpt (int) : Core symmetry (1=full, 4=quarter mirror, 5=quarter rotational)
		"""
		me.rodMap = rodMap
		me.chanMap = chanMap
		me.assemMap = assemMap
		me.symOpt = symOpt

		def _preGenDeckChecks(me):
		"""
		This method is run prior to generating the deck to check for model consistency.

SubKit/build/genDeck.py

+44 −26

Original line number	Diff line number	Diff line
		@@ -91,6 +91,8 @@ def writeDeck(model, filename):
		fluidprops=1
		elif model.fluidProperties=='if97_lookup':
		fluidprops=3
		elif model.fluidProperties=='custom':
		fluidprops=4
		elif model.fluidProperties=='flibe':
		fluidprops=11
		else:
		@@ -660,6 +662,44 @@ def writeDeck(model, filename):
		else:
		raise RuntimeError

		# Data to be written to Card 17
		group17Data=[]
		if model.rodMap:
		group17Data.append("*NGR group number\n")
		group17Data.append(" 17\n")
		group17Data.append("*Card 17.1 - HDF5_NAME VTK_NAME\n")
		group17Data.append(" hdf5_name vtk_name\n")
		group17Data.append("*Card 17.2 - Rod Map Dimensions\n")
		group17Data.append("**TOTRODSROW TOTRODSCOL\n")
		group17Data.append("{0:12d}{1:12d}\n".format(len(model.rodMap), len(model.rodMap[0])))
		group17Data.append("*Card 17.3 - Channel Map Dimensions\n")
		group17Data.append("**TOTCHANSROW TOTCHANSCOL\n")
		group17Data.append("{0:12d}{1:12d}\n".format(len(model.chanMap), len(model.chanMap[0])))
		group17Data.append("*Card 17.4 - Rod Map\n")
		for row in model.rodMap:
		line = ""
		for idx in row:
		line = line+"{:5d}".format(idx)
		line = line+'\n'
		group17Data.append(line)
		group17Data.append("*Card 17.4 - Channel Map\n")
		for row in model.chanMap:
		line = ""
		for idx in row:
		line = line+"{:5d}".format(idx)
		line = line+'\n'
		group17Data.append(line)
		group17Data.append("{assem map}\n")
		group17Data.append("** sym_opt nassem_rows nassem_cols\n")
		group17Data.append("{:10d}{:13d}{:13d}\n".format(model.symOpt, len(model.assemMap), len(model.assemMap[0])))
		group17Data.append("** assembly_map\n")
		for row in model.assemMap:
		line = ""
		for idx in row:
		line = line+"{:5d}".format(idx)
		line = line+'\n'
		group17Data.append(line)

		group19Data = []
		group19Data.append("**NGR\n")
		group19Data.append(" 19\n")
		@@ -701,32 +741,6 @@ def writeDeck(model, filename):
		group19Data.append("** L2AVOID L2APRESS L2ATCOOL L2ATSOLID L2AVL L2AVV L2AVD\n")
		group19Data.append(" 1.000E-04 1.000E-04 1.000E-04 1.000E-04 1.000E-04 1.000E-04 1.000E-04\n")

		# Data to be written to Card 17
		#group17Data=[]
		#group17Data.append("*NGR group number \n")
		#group17Data.append(" 17 \n")
		#group17Data.append("*Card 17.1 - HDF5_NAME VTK_NAME \n")
		#group17Data.append(" hdf5_name vtk_name \n")
		#group17Data.append("*Card 17.2 - Rod Map Dimensions \n")
		#group17Data.append("**TOTRODSROW TOTRODSCOL \n")
		#group17Data.append("{0:d12}{1:d12}\n".format(pinMap.shape(0), pinMap.shape(1)))
		#group17Data.append("*Card 17.3 - Channel Map Dimensions \n")
		#group17Data.append("**TOTCHANSROW TOTCHANSCOL \n")
		#group17Data.append("{0:d12}{1:d12}\n".format(pinMap.shape(0)+1, pinMap.shape(1)+1))
		#group17Data.append("*Card 17.4 - Rod Map \n")
		#for i in range(pinMap.shape(0)):
		# for j in range(pinMap.shape(1)):
		#
		#{assem map}
		#** sym_opt nassem_rows nassem_cols
		# 1 6 6
		#**Assembly Map
		# 1 2 3 4 5 6
		# 7 8 9 10 11 12
		# 13 14 15 16 17 0
		# 18 19 20 21 22 0
		# 23 24 25 26 0 0
		# 27 28 0 0 0 0

		def makeGroupBanner(title):
		assert(isinstance(title, str))
		@@ -783,6 +797,10 @@ def writeDeck(model, filename):
		newDeckLines.append(makeGroupBanner("GROUP 15 - Time Domain Data"))
		for l in group15Data:
		newDeckLines.append(l)
		if group17Data:
		newDeckLines.append(makeGroupBanner("GROUP 17 - Channel/Rod Maps"))
		for l in group17Data:
		newDeckLines.append(l)
		if notrans==1:
		newDeckLines.append(makeGroupBanner("GROUP 19 - convergence Criteria for Steady State Solve"))
		for l in group19Data:

SubKit/process/Hdf5Tools.py

+73 −0

Original line number	Diff line number	Diff line
		@@ -74,6 +74,16 @@ class Hdf5Interface(object):
		nrad = me.getPinNumRadial(pin)
		me.pinInternalMask[0:nsec, 0:nlev, 0:nrad, pin-1] = False

		# Create a mask for surface values from the pin_temp dataset
		if 'STATE_0001/pin_temp' in me.h5:
		me.pinTempSurfMask = np.ones(me.h5['STATE_0001/pin_temp'].shape, dtype=bool)
		me.pinTempSurfMask[:, :, :, :] = True
		for pin in range(1, me.numPin+1):
		nlev = me.getPinNumLev(pin)
		nsec = me.getPinNumSector(pin)
		# Only the surface value is unmasked
		me.pinTempSurfMask[0:nsec, 0:nlev, -1, 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)
		@@ -1022,6 +1032,69 @@ class Hdf5Interface(object):
		data = me.getPinSurfaceDataset('pin_surf_q_tot', state, pin)
		return axial, data[:, sector]

		def getPinWithMaxSurfTemp(me, state=None):
		""" Returns the location with maximum pin surface temperature

		If State is passed, only searches that state point.

		Args:
		state (int) : 1-based state index. Searches all if not passed.

		Returns:
		state (int) : 1-based state index with maximum pin surface temperature.
		pin (int) : 1-based pin index containing maximum surface temperature.
		level (int) : 0-based level index containing maximum pin surface temperature.
		sector (int) : 0-based sector index containing maximum pin surface temperature.
		temp (float) : Surface temperature at max location
		"""
		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_, 'pin_temp')], mask=me.pinTempSurfMask)
		maxval = ma.max()
		if maxVal is None or maxval>maxVal:
		maxVal = maxval
		maxState = state_
		(maxTheta, maxLevel, maxRad, maxPin) = np.unravel_index(ma.argmax(), ma.shape)
		return maxState, maxPin+1, maxLevel, maxTheta, maxVal


		def getPinAxialTempSurf(me, state=None, pin=None, sector=None):
		""" Returns the axial vector mesh and axial surface temp distribution on that mesh for selected pin.

		Args:
		state (int): State number to find pin (1-based). If not provided, picks the state where the
		maximum was experienced.
		pin (int): Pin index to plot (1-based). If not provided, picks the pin where the maximum
		was experienced.
		sector (int): The azimuthal index (0-base). If not provided, picks the maximum value in the
		azimuthal direction.

		Returns:
		axialMesh (list): Axial mesh list
		temp (list): Temperature distribution
		"""
		if state is None: assert(sector is None)
		if state is not None:
		me._stateValid(state)
		if pin is not None: me._pinIndexValid(pin)
		if sector is not None: me._pinSectorValid(pin, sector)

		_state, _pin, _level, _sector, maxVal = me.getPinWithMaxSurfTemp(state)

		if state is not None: _state = state
		if pin is not None: _pin = pin
		if sector is not None: _sector = sector

		data = me.h5['STATE_{:04d}/pin_temp'.format(_state)][()][_sector, 0:me.getPinNumLev(_pin), me.getPinNumRadial(_pin)-1, _pin-1]
		axial = me.h5['CORE/mesh/PIN_{:05d}/axial_nodes'.format(_pin)][()]

		return list(axial), list(data)

		def getChanScalarMesh(me, chan):
		""" Returns the channel scalar mesh including ghost levels

SubKit/process/PlotTools.py

+54 −13

Original line number	Diff line number	Diff line
		@@ -201,23 +201,48 @@ class PlotTools(object):
		else:
		figname_ = figname

		me._plotAxial(list(dnbr), list(axial), 'DNBR [-]', figname_, "DNBR", location, xlim=[0.0, min(dnbr)+10.0])

		fig, ax = plt.subplots()
		plt.title(location)
		ax.grid()
		ax.set_ylabel("Axial [m]")
		ax.set_ylim([axial[0], axial[-1]])
		ax.set_xlabel("DNBR [-]")
		ax.set_xlim([0.0, min(dnbr)+10.0])
		ax.plot(dnbr, axial)
		#ax.legend()
		plt.savefig(figname_)
		plt.close(fig)
		def plotPinAxialTemp(me, state=None, pin=None, theta=None, figname=None):
		""" Makes an axial plot of surface temperature for a pin

		Invalid use case is to pass only theta.

		Args:
		state (int): State number to plot (1-based). If not present, selects state with maximum.
		pin (int): Pin number to plot (1-based). If not present, selects pin with maximum.
		theta (int): Azimuthal index on rod surface to plot (0-based). If not present, selects azimuthal
		location with maximum.
		figname (str): Name of the figure.
		"""
		if ((theta is not None) and ((state is None) and (pin is None))):
		raise RuntimeError("Must pass (state/pin/theta), (state/pin), (state), or None for axial location")

		axial, temp = me.h5obj.getPinAxialTempSurf(state, pin, theta)

		state_, pin_, level_, sector_, maxVal = me.h5obj.getPinWithMaxSurfTemp(state)

		if state is not None:
		state_ = state
		if pin is not None:
		pin_ = pin
		if theta is not None:
		theta_ = theta

		# Make all indices 1-based for I/O
		location = "S{:d}_P{:d}_T{:d}".format(state_, pin_, sector_+1)

		if not figname:
		figname_ = location+"_pin_temp_surf.png"
		else:
		figname_ = figname

		me._plotAxial(temp, axial, 'Temperature [C]', figname_, "Tw", location)

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

		def _plotAxial(me, x, y, xname, figname, labels, title=None, yname="Axial [m]"):
		def _plotAxial(me, x, y, xname, figname, labels, title=None, yname="Axial [m]", xlim=None):
		""" Makes an axial plot

		x (list or dict): x-axis data to plot. Pass a dictionary for multiple datasets.
		@@ -228,6 +253,7 @@ class PlotTools(object):
		a list of equal length to the dicts. It gives the labels for including in the legend.
		title (str): Title for the plot
		yname (str): The name of the y plot
		xlim (list): [min, max] of x-axis in plot.

		"""
		fig, ax = plt.subplots()
		@@ -236,6 +262,10 @@ class PlotTools(object):
		ax.grid()
		ax.set_ylabel(yname)
		ax.set_xlabel(xname)
		if xlim is not None:
		assert(isinstance(xlim, list))
		assert(len(xlim)==2)
		ax.set_xlim(xlim)
		assert(isinstance(y, list) or isinstance(y, dict))
		if isinstance(y, list):
		assert(isinstance(x, list))
		@@ -251,6 +281,7 @@ class PlotTools(object):
		for i, name in enumerate(y.keys()):
		ax.plot(x[name], y[name], label=labels[i])
		addLegend=True
		if addLegend:
		ax.legend()
		plt.savefig(figname)
		plt.close(fig)
		@@ -311,6 +342,16 @@ class PlotTools(object):
		"""
		me._chanPlotMaxWrapper('chan_equilibrium_quality', "Quality [-]", state, ch, figname)

		def plotChAxialTemperature(me, state=None, ch=None, figname=None):
		""" Makes an axial plot of temperature 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_temp', "Temperature [C]", state, ch, figname)

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

SubKit/process/plotChanTemp.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 temperature 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 temperature.")
		parser.add_argument('--ch', type=int, action='append', help="""Channel to plot. Defaults to channel with max temperature.
		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.plotChAxialTemperature(args.state, args.ch, args.figname)

		if __name__=='__main__':
		main()