Fix 3D power normalization. (1cf82871) · Commits · Gurecky, William / SubKit

SubKit/build/CoreBuilder.py

+19 −9

Original line number	Diff line number	Diff line
		@@ -128,14 +128,6 @@ class Core(object):
		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 """
		@@ -152,9 +144,27 @@ class Core(object):
		[z_1, z_2], # axial edit i=1
		[z_2, z_3], ...]
		"""
		z_pairs = np.asarray((me.z_bounds[1:], me.z_bounds[:-1])).T
		z_pairs = np.asarray((me.z_bounds[:-1], me.z_bounds[1:])).T
		return z_pairs

		@property
		def z_power_grid(me):
		"""
		Returns a flattened view of z_bounds_pairs with small adjustments
		made to the cell boundaries. This is for use only when
		specifying a CTF power factor table. The flattened array has
		the form:
		[z_0, z_1-eps, z_1, z_2-eps, z_2, z_3-eps, ...]
		where eps is a small value (eps=1e-6).
		"""
		z_pwr_grid = me.z_bounds_pairs
		# slightly adjust the upper bounds of each axial edit zone
		z_pwr_grid[:, 1] -= 1e-5
		# flatten and return
		z_pwr_grid = z_pwr_grid.flatten()
		z_pwr_grid[-1] = np.cumsum(me.dz)[-1] + 1e-5
		return z_pwr_grid

		# === Delegate to me.model ===
		def __getattr__(me, attr):
		"""

SubKit/build/SquareLatticeLWR_Nodal.py

+34 −14

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


		class NodalChannelGeom(object):
		@@ -1046,7 +1047,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		assert len(me.z_power_grid.shape) == 1
		assert powerDist.shape[0] == me.chMapObj.getAssemDimLen()
		assert powerDist.shape[1] == me.chMapObj.getAssemDimLen()
		assert me.z_power_grid.size == powerDist.shape[2]
		assert me.z_centers.size == powerDist.shape[2]

		# normalize 3d powers
		if wgts is None:
		@@ -1063,12 +1064,22 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		assemID = me.chMapObj.getAssemIndexFromLoc(i, j)
		if assemID:
		powID = "pwr_" + str(assemID)
		asm_axial_power_profile_avg = \
		np.max((1e-12, np.average(powerDist[i, j, :])))
		me.model.addAxialPowerShape(powID, me.z_power_grid, powerDist[i, j, :] / \
		asm_axial_power_profile_avg)
		rods = me.rodsInAssem[assemID]

		asm_axial_power_profile = powerDist[i, j, :]

		# the asm_axial_power_profile is specified at z_centers
		assert me.z_centers.size == asm_axial_power_profile.size
		asm_axial_power_interpolant = \
		interp1d(me.z_centers, asm_axial_power_profile,
		kind='linear', fill_value='extrapolate')

		# interpolate to the z_bounds
		asm_axial_power_ctf = asm_axial_power_interpolant(me.z_power_grid)
		avg_axpwr = np.max((np.trapz(asm_axial_power_ctf, x=me.z_power_grid) / \
		(np.max(me.z_power_grid) - np.min(me.z_power_grid)), 1e-12))

		me.model.addAxialPowerShape(powID, me.z_power_grid, asm_axial_power_ctf / avg_axpwr)
		rods = me.rodsInAssem[assemID]
		for rodID in rods:
		obj_rod = me.model.solidObjects[rodID]
		obj_rod.powID = powID
		@@ -1077,7 +1088,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		# obj_rod_type_id == 2 is guide tube (only present if me.modelGuideTubes=True)
		# obj_rod is type SubKit.build.Solid.HeatedSolid
		# radial power factors are implicit in the 3d power distribution
		obj_rod.power = asm_axial_power_profile_avg
		obj_rod.power = avg_axpwr

		def _setCorePowerShapeNodal3D(me, powerDist, wgts=None, coreStart=0.0):
		"""
		@@ -1106,7 +1117,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[3]
		assert me.z_centers.size == powerDist.shape[3]

		# normalize 3d powers
		if wgts is None:
		@@ -1154,15 +1165,24 @@ 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, :]
		avg_node_axial_power_profile = \
		np.max((1e-12, np.average(node_axial_power_profile)))
		# scale axial power
		me.model.addAxialPowerShape(powID, me.z_power_grid,
		node_axial_power_profile / avg_node_axial_power_profile)

		# the node_axial_power_profile is specified at z_centers
		assert me.z_centers.size == node_axial_power_profile.size
		node_axial_power_interpolant = \
		interp1d(me.z_centers, node_axial_power_profile,
		kind='linear', fill_value='extrapolate')

		# interpolate to the z_bounds
		node_axial_power_ctf = node_axial_power_interpolant(me.z_power_grid)
		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))

		# set axial power profile in the ctf model
		me.model.addAxialPowerShape(powID, me.z_power_grid, node_axial_power_ctf / avg_axpwr)
		obj_rod.powID = powID
		obj_rod_type_id = me.model.solidTypeIDs[rodID]
		# obj_rod_type_id == 1 is heated rod
		# obj_rod_type_id == 2 is guide tube (only present if me.modelGuideTubes=True)
		# obj_rod is type SubKit.build.Solid.HeatedSolid
		# radial power factors are implicit in the 3d power distribution
		obj_rod.power = avg_node_axial_power_profile
		obj_rod.power = avg_axpwr

setup.py

+1 −0

Original line number	Diff line number	Diff line
		@@ -30,6 +30,7 @@ setup(name='SubKit',
		'mock',
		'matplotlib<=2.2.4',
		'sympy',
		'scipy',
		'numpy<=1.16.4',
		'h5py'
		]

tests/nodal_p7_pow3d/deck.inp.gold

+79419 −41591

File changed.

Preview size limit exceeded, changes collapsed.