Merge remote-tracking branch 'willfork/set_3dpower_dist' into calc_nodal_powers (420dc34f) · Commits · Phathanapirom, Birdy / SubKit

SubKit/build/CoreBuilder.py

+48 −1

Original line number	Diff line number	Diff line
		@@ -59,6 +59,53 @@ class Core(object):
		""" Returns a list of the channel IDs of the channels in the core region."""
		return me.core_ch

		### Core Attributes ###
		@property
		def dz(me):
		""" Mesh partition heights """
		return me._dz

		@dz.setter
		def dz(me, dz_in):
		"""
		Set mesh partition heights.

		Args:
		dz_in (list of float): A list of the heights of each axial level in the core.
		"""
		me._dz = np.asarray(dz_in)
		assert len(me._dz.shape) == 1
		assert np.min(me._dz) >= 0.0
		assert np.max(me._dz) <= 1.0e10

		@property
		def z_centers(me):
		""" Axial level z centers. Has length = len(z_bounds) - 1 """
		z_center_grid = np.average(me.z_bounds_pairs, axis=1)
		return z_center_grid

		@property
		def z_power_grid(me):
		""" Adjust z center endpoints to span entire axial core height """
		z_pwr_grid = me.z_centers
		z_pwr_grid[0] = 0.0
		z_pwr_grid[-1] = np.cumsum(me.dz)[-1] + 0.0001
		return z_pwr_grid

		@property
		def z_bounds(me):
		""" Axial level z boundaries """
		z_bound = np.cumsum(me.dz)
		z_bound = np.insert(z_bound, 0, 0.0)
		return z_bound

		@property
		def z_bounds_pairs(me):
		""" Axial level z boundary pairs ((zstart_i, zend_i), ...) """
		z_pairs = np.asarray((me.z_bounds[1:], me.z_bounds[:-1])).T
		return z_pairs


		class MSRE(Core):
		""" Core extension for building core geometry map for the MSRE.

		@@ -73,7 +120,7 @@ class MSRE(Core):
		a blank spot in the map. Must be square in shape. Provide the full symmetry map.
		sym_opt (int): Enter either 1 for full core geometry or 4 for quarter symmetry geometry.
		grBlockWidth (float): The width of a graphite block.
		dz (float): A list of the heights of each axial level in the core.
		dz (list of float): A list of the heights of each axial level in the core.
		startChanIndex (int): The channels must have an integer index when using this map builder. This is the
		index of the first channel in the core region.
		chLen (float): The channels in the MSRE are football shaped. This is the length of that shape [m].

SubKit/build/Model.py

+100 −4

Original line number	Diff line number	Diff line
		@@ -31,7 +31,7 @@ class Model:
		me.cha = {}
		me.chb = {}
		# List of all gap objets in the model
		me.gaps = []
		me.gaps = {}
		# Takes the solid ID (user defined) and returns the solid type ID (user defined)
		me.solidTypeIDs = {}
		# Takes the solid type ID (user defined) and returns the solid geometry object it refers to
		@@ -136,6 +136,9 @@ class Model:
		me.setInitialConditionsCalled = False
		me.setInitialConditionsMassFluxCalled = False

		# Single-phase turbulent mixing coefficient
		me.beta = 0.037


		def setTitle(me, title):
		"""Set the title of the model.
		@@ -146,6 +149,16 @@ class Model:
		"""
		me.title = str(title)

		def setModelOptions(me, beta=None):
		""" Set various modeling options

		Args:
		beta (float): Single-phase turbulent mixing coefficietn
		"""

		if beta is not None:
		me.beta = beta

		def setFluidProperties(me, fluidprops):
		"""Set the fluid property table to use.

		@@ -293,7 +306,7 @@ class Model:
		else:
		me.chb[chID].append(chb)

		def setGap(me, ii, jj, width, length, k=DEFAULT_GAP_FORM_LOSS, mult=1.0):
		def setGap(me, ii, jj, width, length, k=DEFAULT_GAP_FORM_LOSS, mult=1.0, gapID=None, gapBelow=0, gapAbove=0):
		""" Defines a lateral connection between two channels in the same axial section

		Args:
		@@ -303,9 +316,25 @@ class Model:
		length (float): Lateral distance between the centroids of the two connected channels
		k (float): Form loss applied to the lateral momentum equation solved in the gap
		mult (float): Number of physical gaps in the model represented by this gap
		gapID (int) : ID of the gap. If not assigned, an unused one will be created and returned.
		gapBelow (int): ID of gap below this one. This is used for transporting lateral momentum
		accross section boundaries and therefore must be specified for multisection models where
		gaps are connected. Enter 0 or omit if no gap below.
		gapAbove (int): ID of gap above this one. Enter 0 or omit if no gap above.

		"""
		me.gaps.append(Gap(ii, jj, width, length, k, mult))
		if gapID is None:
		gapIDs = sorted(me.gaps.keys())
		if len(gapIDs)>0:
		gapID = gapIDs[-1]+1
		else:
		gapID = 1
		else:
		assert isinstance(gapID, int)

		me.gaps[gapID] = Gap(ii, jj, width, length, k, mult, gapBelow, gapAbove)

		return gapID

		def solidQuickAdd(me, solidID, geoObj, solidObj=None):
		""" Used to add a single instance of a solid object
		@@ -945,6 +974,22 @@ class Model:
		me.channelsInDomain[domain].append(IDsInDomain)
		me.solver = 5

		def getSectionOfGap(me, gapID):
		""" Returns the secID of the section that owns the passed gapID

		Args:
		gapID (int) : Gap ID

		"""
		assert gapID in me.gaps
		sec1 = me.getSectionOfChannel(me.gaps[gapID].ii)
		sec2 = me.getSectionOfChannel(me.gaps[gapID].jj)
		if sec1!=sec2:
		raise RuntimeError("Gap {:d} connects Channels {:d} and {:d} which are in different sections "+
		"{:d} and {:d}. This is not allowed".format(gapID, me.gaps[gapID].ii, me.gaps[gapID].jj, sec1, sec2))
		return sec1


		def getSectionOfChannel(me, chID):
		""" Returns the secID of the section that owns the passed chID

		@@ -1246,6 +1291,8 @@ class Model:

		me._rodConnectionsValid()

		me._gapChecks()



		def _rodsInValidSections(me):
		@@ -1352,6 +1399,51 @@ class Model:
		if chID not in me.cha[below]:
		me.cha[below].append(chID)

		def _gapChecks(me):
		""" Ensure gaps defined in model are valid"""
		if len(me.gaps)==0:
		return
		# Gap IDs must start at 1 and run sequentially to num gaps
		gapIDs = sorted(me.gaps.keys())
		errMsg = "Gap IDs must start at 1 and run continuously to total number of gaps"
		if gapIDs[0]!=1:
		raise RuntimeError(errMsg)
		for gapCount in range(1, len(gapIDs)):
		if gapIDs[gapCount]!=gapIDs[gapCount-1]+1:
		raise RuntimeError(errMsg)
		# gaps that connect across section bounds must be connected to valid indices
		for gapID in me.gaps.keys():
		gapSection = me.getSectionOfGap(gapID)
		gapBelow = me.gaps[gapID].gapBelow
		if gapBelow!=0:
		gapSectionBelow = me.getSectionOfGap(gapBelow)
		# Make sure the section of the gap below is actually the section below this gap's section
		if gapSectionBelow!=gapSection-1:
		raise RuntimeError("Gap {:d} connects to Gap {:d} below, but Gap {:d} is in Section {:d} and "+
		"Gap {:d} is in Section {:d}, which are not correctly axially aligned".format(gapID,
		gapBelow, gapID, gapSection, gapBelow, gapSectionBelow))
		# If the gap below does not have gapAbove set to be consistent with this one, set it
		if me.gaps[gapBelow].gapAbove==0:
		me.gaps[gapBelow].gapAbove = gapID
		# If the gap below was set, but is not equal to this gapID, throw an error
		if me.gaps[gapBelow].gapAbove!=gapID:
		raise RuntimeError("Gap {:d} had gapAbove set to {:d}, but this is not consistent with what "+
		"gapBelow was set to for {:d}".format(gapBelow, me.gaps[gapBelow].gapAbove, gapID))
		gapAbove = me.gaps[gapID].gapAbove
		if gapAbove!=0:
		gapSectionAbove = me.getSectionOfGap(gapAbove)
		if gapSectionAbove!=gapSection+1:
		raise RuntimeError("Gap {:d} connects to Gap {:d} above, but Gap {:d} is in Section {:d} and "+
		"Gap {:d} is in Section {:d}, which are not correctly axially aligned".format(gapID,
		gapAbove, gapID, gapSection, gapAbove, gapSectionAbove))
		# If the gap above does not have gapAbove set to be consistent with this one, set it
		if me.gaps[gapAbove].gapBelow==0:
		me.gaps[gapAbove].gapBelow = gapID
		# If the gap above was set, but is not equal to this gapID, throw an error
		if me.gaps[gapAbove].gapBelow!=gapID:
		raise RuntimeError("Gap {:d} had gapBelow set to {:d}, but this is not consistent with what "+
		"gapAbove was set to for {:d}".format(gapAbove, me.gaps[gapAbove].gapBelow, gapID))

		def _validateBoundaryConditions(me):
		""" Ensure boundary conditions were setup correctly"""
		# Ensure all were defined where expected
		@@ -1409,12 +1501,16 @@ class Gap(object):
		length (float): Lateral distance between the centroids of the two connected channels
		k (float): Form loss applied to the lateral momentum equation solved in the gap
		mult (float): Number of physical gaps in the model represented by this gap
		gapBelow (int) : ID of gap below this one. Enter 0 or omit if no gap below.
		gapAbove (int) : ID of gap above this one. Enter 0 or omit if no gap above.

		"""
		def __init__(me, ii, jj, width, length, k=DEFAULT_GAP_FORM_LOSS, mult=1.0):
		def __init__(me, ii, jj, width, length, k=DEFAULT_GAP_FORM_LOSS, mult=1.0, gapBelow=0, gapAbove=0):
		me.ii = ii
		me.jj = jj
		me.width = width
		me.length = length
		me.k = k
		me.mult = mult
		me.gapBelow = gapBelow
		me.gapAbove = gapAbove

SubKit/build/SquareLatticeLWR_Nodal.py

+39 −36

Original line number	Diff line number	Diff line
		@@ -834,15 +834,15 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		""" Use to set a power shape in the core region.

		Args:
		radialPower (func): A function that takes non-dimensional radial location and returns
		a nondimensional power factor. If not provided, a uniform radial power distribution
		radialPower (np_2darray or list): A two dimensional list or numpy array that contains radial
		nondimensional power factors. If not provided, a uniform radial power distribution
		is assumed.
		axialPower (np_1darray or list): A one dimensional list or numpy array that contains
		nondimensional axial multiplier on power. If not provided, a uniform axial power distribution
		is assumed.
		axialPower (func): A function that takes axial location and returns a nondimensional
		multiplier on power. If not provided, a uniform axial power distribution is assumed.
		coreStart (float): The axial location of the start of the core. In CTF, the bottom of the
		model is zero. If the core does not start at the bottom of the model, this must be provided
		so the power tables are correctly defined in CTF. [m]

		"""
		# Go back through the pins and assign radial power factors to all
		# Just define one axial shape map
		@@ -896,34 +896,29 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		for i, obj in enumerate(me.model.solidObjects.values()):
		obj.power = radP[i]

		@property
		def z_centers(me):
		""" Axial level z centers. Has length = len(z_bounds) - 1 """
		z_center_grid = np.average(me.z_bounds_pairs, axis=1)
		return z_center_grid

		@property
		def z_power_grid(me):
		""" Adjust z center endpoints to span entire axial core height """
		z_pwr_grid = me.z_centers
		z_pwr_grid[0] = 0.0
		z_pwr_grid[-1] = np.cumsum(me.dz)[-1] + 0.0001
		return z_pwr_grid

		@property
		def z_bounds(me):
		""" Axial level z boundaries """
		z_bound = np.cumsum(me.dz)
		z_bound = np.insert(z_bound, 0, 0.0)
		return z_bound

		@property
		def z_bounds_pairs(me):
		""" Axial level z boundary pairs ((zstart_i, zend_i), ...) """
		z_pairs = np.asarray((me.z_bounds[1:], me.z_bounds[:-1])).T
		return z_pairs

		def setCorePowerShape3D(me, powerDist, wgts=None, coreStart=0.0):
		"""
		Set a 3D power distribution for the core.
		Args:
		powerDist (numpy.ndarray):
		3 dim array. Contains power distribution with indexing
		(asm_i, asm_j, asm_axial)
		or
		4 dim array. Contains power distribution with indexing
		(asm_i, asm_j, node_id, asm_axial)
		wgts (numpy.ndarray): Same shape as powerDist. Weights of powers (optional)
		coreStart (float): The axial location of the start of the core. In CTF, the bottom of the
		model is zero. If the core does not start at the bottom of the model, this must be provided
		so the power tables are correctly defined in CTF. [m]
		"""
		if len(powerDist.shape) == 3:
		me._setCorePowerShapeAsm3D(powerDist, wgts, coreStart)
		elif len(powerDist.shape) == 4:
		me._setCorePowerShapeNodal3D(powerDist, wgts, coreStart)
		else:
		raise RuntimeError("Invalid power distribution shape: %s" % str(powerDist.shape))

		def _setCorePowerShapeAsm3D(me, powerDist, wgts=None, coreStart=0.0):
		"""
		Set a 3D power distribution for the core. Each assembly can have it's own
		axial power shape.
		@@ -971,7 +966,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		# radial power factors are implicit in the 3d power distribution
		obj_rod.power = 1.0

		def setCorePowerShapeNodal3D(me, powerDist, wgts=None, coreStart=0.0):
		def _setCorePowerShapeNodal3D(me, powerDist, wgts=None, coreStart=0.0):
		"""
		Set a 3D power distribution for the core. Each assembly can have it's own
		axial power shape.
		@@ -979,6 +974,14 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		Args:
		powerDist (numpy.ndarray): 4 dim array. Contains power distribution with indexing
		(asm_i, asm_j, node_id, asm_axial)
		node_id is indexed from left to right, top to bottom:
		inter-asm node_id layout:
		---------
		\| 0 \| 1 \|
		---------
		\| 2 \| 3 \|
		---------

		wgts (numpy.ndarray): Same shape as powerDist. Weights of powers (optional)
		coreStart (float): The axial location of the start of the core. In CTF, the bottom of the
		model is zero. If the core does not start at the bottom of the model, this must be provided
		@@ -990,7 +993,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		assert powerDist.shape[0] == me.chMapObj.getAssemDimLen()
		assert powerDist.shape[1] == me.chMapObj.getAssemDimLen()
		assert powerDist.shape[2] == 4
		assert me.z_power_grid.size == powerDist.shape[2]
		assert me.z_power_grid.size == powerDist.shape[3]

		# normalize 3d powers
		if wgts is None:
		@@ -1022,10 +1025,10 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		obj_rod_type_id = me.model.solidTypeIDs[rodID]
		if obj_rod.x > centroid_xy[0] and obj_rod.y > centroid_xy[1]:
		# north east
		nodeID = 0
		nodeID = 1
		elif obj_rod.x < centroid_xy[0] and obj_rod.y > centroid_xy[1]:
		# north west
		nodeID = 1
		nodeID = 0
		elif obj_rod.x < centroid_xy[0] and obj_rod.y < centroid_xy[1]:
		# south west
		nodeID = 2

SubKit/build/genDeck.py

+6 −5

Original line number	Diff line number	Diff line
		@@ -156,12 +156,13 @@ def writeDeck(model, filename):
		group3Data.append("** K IK JK GAP LNGT WKR FWAL IGPB IGPA FACT IGAP JGAP IGAP JGAP IGAP JGAP\n")
		group3Data.append("*Card 3.3\n")
		group3Data.append("**GMULT ETNR\n")
		for k, gap in enumerate(model.gaps):
		group3Data.append("{0:6d}{1:6d}{2:6d}{3:16.6e}{4:16.6e}{5:15.5e} 0.0 0 0 1.0 0 0 0 0 0 0\n".format(k+1, chMap.getIndex(gap.ii), chMap.getIndex(gap.jj), gap.width, gap.length, gap.k))
		for k in model.gaps.keys():
		gap = model.gaps[k]
		group3Data.append("{:6d}{:6d}{:6d}{:16.6e}{:16.6e}{:15.5e} 0.0{:5d}{:5d} 1.0 0 0 0 0 0 0\n".format(k, chMap.getIndex(gap.ii), chMap.getIndex(gap.jj), gap.width, gap.length, gap.k, gap.gapBelow, gap.gapAbove))
		group3Data.append("{:6.4f} 0.0\n".format(gap.mult))
		group3Data.append("*Card 3.3.5\n")
		for k, gap in enumerate(model.gaps):
		group3Data.append("{0:6d}{1:>6s}\n".format(k+1, 'x'))
		for k in model.gaps.keys():
		group3Data.append("{0:6d}{1:>6s}\n".format(k, 'x'))
		group3Data.append("*Card 3.4\n")
		group3Data.append("**NLGP\n")
		group3Data.append(" 0\n")
		@@ -629,7 +630,7 @@ def writeDeck(model, filename):
		group12Data.append(" 12\n")
		group12Data.append("*Card 12.2\n")
		group12Data.append("** AAAK BETA THETM\n")
		group12Data.append(" 0.140E+01 0.500E-02 0.500E+01\n")
		group12Data.append(" 0.140E+01{:10.3e} 0.500E+01\n".format(model.beta))


		group13Data.append("**Card 13.4\n")

SubKit/process/PlotTools.py

+29 −41

Original line number	Diff line number	Diff line
		@@ -322,48 +322,32 @@ class PlotTools(object):
		labels.append("ch {:4d}".format(c))
		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
		def plotChAxial(me, field, state=None, ch=None, figname=None):
		""" Makes an axial plot of channel data

		Args:
		field (str): The name of the data to plot. Select from the following---'massflux' for total,
		massflux, 'temp' for liquid temperature, 'eq_quality' for equilibrium quality, 'drop_velocity'
		for droplet axial velocity, 'vap_velocity' for vapor axial velocity, 'liq_velocity' for liquid
		axial velocity, or 'void' for vapor void.
		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 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

		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)

		chData = {'massflux': ['massflux_tot', 'Mass flux [kg/m**2/s]'],
		'temp': ['chan_temp', 'Temperature [C]'],
		'eq_quality': ['chan_equilibrium_quality', 'Quality [-]'],
		'drop_velocity': ['chan_velocity_drp', 'Velocity [m/s]'],
		'vap_velocity': ['chan_velocity_vap', 'Velocity [m/s]'],
		'liq_velocity': ['chan_velocity_liq', 'Velocity [m/s]'],
		'void': ['chan_void_vap', 'Void [-]']
		}
		if field in chData:
		me._chanPlotMaxWrapper(chData[field][0], chData[field][1], state, ch, figname)
		else:
		raise RuntimeError("Selected dataset: {:s} is not a valid field name".format(field))

		def plotPinCoupling(me, pin, dataset, state=None, figname=None, axialBounds=None, maxTemp=None):
		def plotPinCoupling(me, pin, dataset, state=None, figname=None, axialBounds=None, minData=None, maxData=None):
		""" Makes a contour plot of rod surface data

		Passed dataset must be a coupling dataset (has prefix of "coupling")
		@@ -376,7 +360,8 @@ class PlotTools(object):
		axialBounds (float) : List of 2 values. Specify min and max axial locations to include
		in plot. Specify "None" to use the dataset min or max. First value is the min and
		second is the max.
		maxTemp (float) : Maximum temp to show in contour plot [C]
		minData (float) : Minimum value to show in contour plot
		maxData (float) : Maximum value to show in contour plot

		"""
		if state is None:
		@@ -405,11 +390,14 @@ class PlotTools(object):
		if axialBounds_[1] is None:
		axialBounds_[1] = max(axial)

		if maxTemp is None:
		maxTemp_ = np.max(data)
		if maxData is None:
		maxData_ = np.max(data)
		else:
		maxData_ = maxData
		if minData is None:
		minData_ = np.min(data)
		else:
		maxTemp_ = maxTemp
		minTemp_ = np.min(data)
		minData_ = minData

		fig, ax = plt.subplots()
		plt.title(location)
		@@ -419,7 +407,7 @@ class PlotTools(object):
		#v = np.linspace(minTemp_, maxTemp_, 7, endpoint=True)
		#c1 = ax.contourf(azimuthal, axial, data, v)
		#cbar = plt.colorbar(c1)
		im = ax.pcolor(azimuthal, axial, data, vmin=minTemp_, vmax=maxTemp_)
		im = ax.pcolor(azimuthal, axial, data, vmin=minData_, vmax=maxData_)
		ax.grid()
		cbar = fig.colorbar(im)
		dataset = dataset.replace('_',' ')