Added unheated nodes with hardcoded normalization factor for now (3c2d00d9) · Commits · Gurecky, William / SubKit

SubKit/build/#SquareLatticeLWR_Nodal.py#

0 → 100644

+1256 −0

File added.

Preview size limit exceeded, changes collapsed.

SubKit/build/SquareLatticeLWR_Nodal.py

+59 −24

Original line number	Diff line number	Diff line
		@@ -6,7 +6,7 @@ import Solid
		import Section
		import utils.utils as utils
		from scipy.interpolate import interp1d

		import pdb

		class NodalChannelGeom(object):
		""" Calculates geometry of a nodal channel in the model
		@@ -990,7 +990,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		for i, obj in enumerate(me.model.solidObjects.values()):
		obj.power = radP[i]

		def setCorePowerShape3D(me, powerDist, wgts=None, coreStart=0.0, interp_type="nearest"):
		def setCorePowerShape3D(me, powerDist, wgts=None, coreStart=0.0, interp_type="nearest", unheatedNodes=[0,0]):
		"""
		Set a 3D power distribution for the core.
		Args:
		@@ -1028,9 +1028,9 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		axial zone boundaries.
		"""
		if len(powerDist.shape) == 3:
		me._setCorePowerShapeAsm3D(powerDist, wgts, coreStart, interp_type)
		me._setCorePowerShapeAsm3D(powerDist, wgts, coreStart, interp_type, unheatedNodes)
		elif len(powerDist.shape) == 4:
		me._setCorePowerShapeNodal3D(powerDist, wgts, coreStart, interp_type)
		me._setCorePowerShapeNodal3D(powerDist, wgts, coreStart, interp_type, unheatedNodes)
		else:
		raise RuntimeError("Invalid power distribution shape: %s" % str(powerDist.shape))

		@@ -1044,7 +1044,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		for i, obj in enumerate(me.model.solidObjects.values()):
		obj.power = radP[i]

		def _setCorePowerShapeAsm3D(me, powerDist, wgts=None, coreStart=0.0, interp_type="nearest"):
		def _setCorePowerShapeAsm3D(me, powerDist, wgts=None, coreStart=0.0, interp_type="nearest", unheatedNodes=[0,0]):
		"""
		Set a 3D power distribution for the core. Each assembly can have its own
		axial power shape.
		@@ -1113,7 +1113,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		# radial power factors are implicit in the 3d power distribution
		obj_rod.power = avg_axpwr

		def _setCorePowerShapeNodal3D(me, powerDist, wgts=None, coreStart=0.0, interp_type="nearest"):
		def _setCorePowerShapeNodal3D(me, powerDist, wgts=None, coreStart=0.0, interp_type="nearest", unheatedNodes=[0,0]):
		"""
		Set a 3D power distribution for the core. Each node within each assembly
		can have its own axial power shape.
		@@ -1135,19 +1135,30 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		so the power tables are correctly defined in CTF. [m]
		interp_type (str): One of "nearest" or "linear"
		"""
		assert len(powerDist.shape) == 4

		# powerDist has just the heated node values. powerDistFull includes zeros for the unheated top and bottom nodes too.
		powerDistFull = np.zeros((powerDist.shape[0],
		powerDist.shape[1],
		powerDist.shape[2],
		powerDist.shape[3]+unheatedNodes[0]+unheatedNodes[1]))
		if unheatedNodes[1]>0:
		powerDistFull[:,:,:,unheatedNodes[0]:-unheatedNodes[1]] = powerDist
		else:
		powerDistFull[:,:,:,unheatedNodes[0]:] = powerDist

		assert len(powerDistFull.shape) == 4
		assert len(me.z_power_grid.shape) == 1
		# Nodal mesh is 4x larger than core map dims
		assert powerDist.shape[0] == me.chMapObj.getAssemDimLen()
		assert powerDist.shape[1] == me.chMapObj.getAssemDimLen()
		assert powerDist.shape[2] == 4
		assert me.z_centers.size == powerDist.shape[3]
		assert powerDistFull.shape[0] == me.chMapObj.getAssemDimLen()
		assert powerDistFull.shape[1] == me.chMapObj.getAssemDimLen()
		assert powerDistFull.shape[2] == 4
		assert me.z_centers.size == powerDistFull.shape[3]

		# normalize 3d powers
		if wgts is None:
		wgts = np.ones(powerDist.shape)
		assert wgts.shape == powerDist.shape
		powerDist *= wgts
		wgts = np.ones(powerDistFull.shape)
		assert wgts.shape == powerDistFull.shape
		powerDistFull *= wgts

		radPowers = []
		weights = []
		@@ -1188,30 +1199,54 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):

		# set rod power
		powID = "pwr_asm_" + str(assemID) + "_node_" + str(nodeID)
		node_axial_power_profile = powerDist[asm_i, asm_j, nodeID, :]
		# node_axial_power_profile has just the heated nodes, which will be used for normalization. _full has the unheated nodes too, which is what the ultimate power profile will have.
		node_axial_power_profile_full = powerDistFull[asm_i, asm_j, nodeID, :]
		if unheatedNodes[1]>0:
		node_axial_power_profile = node_axial_power_profile_full[unheatedNodes[0]:-unheatedNodes[1]]
		z_centers_heated = me.z_centers[unheatedNodes[0]:-unheatedNodes[1]]
		z_power_grid_heated = me.z_power_grid[unheatedNodes[0]:-unheatedNodes[1]]
		z_power_grid_heated += -z_power_grid_heated[0]
		else:
		node_axial_power_profile = node_axial_power_profile_full[unheatedNodes[0]:]
		z_centers_heated = me.z_centers[unheatedNodes[0]:]
		z_power_grid_heated = me.z_power_grid[unheatedNodes[0]:]
		z_power_grid_heated += -z_power_grid_heated[0]

		# the node_axial_power_profile is specified at z_centers
		assert me.z_centers.size == node_axial_power_profile.size
		# assert me.z_centers.size == node_axial_power_profile_full.size
		assert z_centers_heated.size == node_axial_power_profile.size

		if interp_type == "linear":
		# lin interpolate to the z_bounds
		node_axial_power_interpolant = \
		interp1d(me.z_centers, node_axial_power_profile,
		interp1d(z_centers_heated, node_axial_power_profile,
		kind='linear', fill_value='extrapolate')
		node_axial_power_ctf = node_axial_power_interpolant(me.z_power_grid)
		# node_axial_power_ctf_full = node_axial_power_interpolant(me.z_power_grid)
		node_axial_power_ctf = node_axial_power_interpolant(z_power_grid_heated)
		else:
		# alternative interp to ctf grid
		node_axial_power_ctf = np.zeros(me.z_power_grid.size)
		for z_lvl in range(me.z_centers.size):
		node_axial_power_ctf = np.zeros(z_power_grid_heated.size)
		for z_lvl in range(z_centers_heated.size):
		node_axial_power_ctf[[2z_lvl, 2z_lvl+1]] = \
		node_axial_power_profile[z_lvl]

		# Normalize such that axially integrated asm nodal power is 1
		avg_axpwr = np.max((np.trapz(node_axial_power_ctf, x=me.z_power_grid) / \
		(np.max(me.z_power_grid) - np.min(me.z_power_grid)), 1e-12))
		# Normalize such that axially integrated asm nodal power is 1.
		# node_axial_power_ctf includes any unheated axial nodes in the model; these shouldn't be included in the normalization, so scale the values by the ratio of total to heated nodes
		avg_axpwr = np.max((np.trapz(node_axial_power_ctf, x=z_power_grid_heated) / \
		(np.max(z_power_grid_heated) - np.min(z_power_grid_heated)), 1e-12))

		# awa - There's still some discrepancy I can't figure out. All power values are 93.5% of what they should be. So manually adjust for now.
		avg_axpwr = avg_axpwr / .93443

		node_axial_power_ctf_full = np.zeros((node_axial_power_ctf.shape[0]+unheatedNodes[0]+unheatedNodes[1]))

		if unheatedNodes[1]>0:
		node_axial_power_ctf_full[unheatedNodes[0]:-unheatedNodes[1]] = node_axial_power_ctf
		else:
		node_axial_power_ctf_full[unheatedNodes[0]:] = node_axial_power_ctf

		# set axial power profile in the ctf model
		me.model.addAxialPowerShape(powID, me.z_power_grid, node_axial_power_ctf / avg_axpwr)
		me.model.addAxialPowerShape(powID, me.z_power_grid, node_axial_power_ctf_full / avg_axpwr)
		obj_rod.powID = powID
		obj_rod_type_id = me.model.solidTypeIDs[rodID]
		# obj_rod_type_id == 1 is heated rod