move axial attributes to base class (068a023e) · Commits · Phathanapirom, Birdy / SubKit

SubKit/build/CoreBuilder.py

+40 −0

Original line number	Diff line number	Diff line
		@@ -59,6 +59,46 @@ 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):
		return me._dz

		@dz.setter
		def dz(me, dz_in):
		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.

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

tests/nodal_p7_pow3d/buildDeck.py

+32 −18

Original line number	Diff line number	Diff line
		@@ -14,7 +14,7 @@ import sys
		sys.path.insert(0, '../../')


		def main(power_3d_flag=True):
		def main(power_3d_flag=1):

		dFuelOuter = 0.475e-2*2 # m
		dFuelInner = 0.418e-2*2 # m
		@@ -199,7 +199,8 @@ def main(power_3d_flag=True):
		2: cdb.Tube(d_outer=dGuideTubeOuter, d_inner=dGuideTubeInner)} # Guide tube

		# Radial power shape
		radPow = [[0.00, 0.00, 0.00, 0.00, 0.75, 0.76, 0.77, 0.78, 0.77, 0.76, 0.75, 0.00, 0.00, 0.00, 0.00],
		radPow = [[0.00, 0.00, 0.00, 0.00, 0.75, 0.76, 0.77, 0.78,
		0.77, 0.76, 0.75, 0.00, 0.00, 0.00, 0.00],
		[0.00, 0.00, 0.00, 0.75, 1.10, 0.93, 1.10, 0.91,
		1.10, 0.93, 1.10, 0.75, 0.00, 0.00, 0.00],
		[0.00, 0.00, 0.74, 1.12, 0.93, 1.17, 0.97, 1.18,
		@@ -226,7 +227,8 @@ def main(power_3d_flag=True):
		0.97, 1.17, 0.93, 1.12, 0.74, 0.00, 0.00],
		[0.00, 0.00, 0.00, 0.75, 1.10, 0.93, 1.10, 0.91,
		1.10, 0.93, 1.10, 0.75, 0.00, 0.00, 0.00],
		[0.00, 0.00, 0.00, 0.00, 0.75, 0.76, 0.77, 0.78, 0.77, 0.76, 0.75, 0.00, 0.00, 0.00, 0.00]]
		[0.00, 0.00, 0.00, 0.00, 0.75, 0.76, 0.77,
		0.78, 0.77, 0.76, 0.75, 0.00, 0.00, 0.00, 0.00]]

		builder = SquareLatticeLWR_Nodal(model, coreMap, pwr17x17, pitch, dz, solidGeoms, assemPitch=assemPitch,
		baffleGap=baffleGap, assemLosses=formLoss, assemLossLevels=gridLevels, rodDomains=rodDomains,
		@@ -256,28 +258,40 @@ def main(power_3d_flag=True):

		coreHeight = np.sum(np.array(dz))

		def axPow(z):
		return np.cos((z-coreHeight/2)*np.pi/2/coreHeight)
		# Power as function of axial position
		axPow = lambda z: np.cos((z-coreHeight/2)*np.pi/2/coreHeight)

		if power_3d_flag:
		# set 3d power distribution
		if power_3d_flag == 1:
		radPow = np.asarray(radPow)
		axialMesh = builder.z_centers
		radPow_arr = np.array(radPow)
		radPow_arr_full = np.ones(
		(radPow_arr.shape[0], radPow_arr.shape[1], len(axialMesh)))
		pow3D = np.ones(
		(radPow.shape[0], radPow.shape[1], len(axialMesh)))
		for z_i in range(len(axialMesh)):
		radPow_arr_full[:, :, z_i] = radPow_arr_full[:,
		:, z_i] * radPow_arr
		axPow_arr = axPow(np.array(axialMesh))
		powerDist3d = radPow_arr_full * axPow_arr
		pow3D[:, :, z_i] = pow3D[:, :, z_i] * radPow
		axPow = axPow(np.array(axialMesh))
		powerDist3d = pow3D * axPow
		builder.setCorePowerShape3D(powerDist3d)
		model.generateModel('deck.inp')
		model.generateModel('asm3dpow_deck.inp')
		elif power_3d_flag == 2:
		radPow = np.asarray(radPow)
		axialMesh = builder.z_centers
		pow3D = np.ones(
		(radPow.shape[0], radPow.shape[1], 4, len(axialMesh)))
		for z_i in range(len(axialMesh)):
		for node_j in range(4):
		pow3D[:, :, node_j, z_i] = pow3D[:, :, node_j, z_i] * radPow
		axPow = axPow(np.array(axialMesh))
		powerDist3d = pow3D * axPow
		builder.setCorePowerShape3D(powerDist3d)
		model.generateModel('nodal3dpow_deck.inp')
		else:
		builder.setCorePowerShape(radialPower=radPow, axialPower=axPow)
		model.generateModel('orig_deck.inp')


		if __name__ == '__main__':
		# build 2 decks which should produce nearly the same CTF TH solution output
		main(power_3d_flag=True)
		main(power_3d_flag=False)
		# build 3 decks which should produce nearly the same CTF TH solution output
		main(power_3d_flag=0)
		main(power_3d_flag=1)
		main(power_3d_flag=2)