Allow user to set discrete 3D power distribution for nodal models (6ac43004) · Commits · Gurecky, William / SubKit

SubKit/build/CoreBuilder.py

+125 −10

Original line number	Diff line number	Diff line
		@@ -21,8 +21,11 @@ import Section

		class Core(object):
		""" Base class for building core maps."""
		def __init__(me):
		pass
		def __init__(me, model, dz, core_sym):
		me.model = model
		me.dz = dz
		me.core_sym = core_sym
		me._isinit = True

		def getSection(me):
		""" Returns the Section object that contains the core. """
		@@ -40,6 +43,10 @@ class Core(object):
		""" See overriding procedure doc"""
		raise NotImplementedError()

		def setEditOptions(me, args, *kwrags):
		""" Set output edit options """
		raise NotImplementedError()

		def getCoreMaxRadialBound(me):
		""" Returns the furthest point in the model from the core origin."""
		coreRadius = 0.0
		@@ -59,6 +66,117 @@ class Core(object):
		""" Returns a list of the channel IDs of the channels in the core region."""
		return me.core_ch

		### Core Base Attributes ###
		@property
		def isinit(me):
		""" Initialization boolian flag. """
		return me._isinit

		@property
		def model(me):
		""" Model container """
		return me._model

		@model.setter
		def model(me, my_model):
		"""
		Set Model container
		Args:
		my_model: instance of Model.Model
		"""
		assert isinstance(my_model, Model.Model)
		me._model = my_model

		@property
		def core_sym(me):
		"""
		Core symmetry property.
		1 == full sym
		4 == qtr sym
		"""
		return me._core_sym

		@core_sym.setter
		def core_sym(me, core_sym_in):
		"""
		Set the core symmetry.
		"""
		assert core_sym_in in [1, 4]
		me._core_sym = core_sym_in

		@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_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 edit z boundary pairs. A 2D array consisting of
		location of each edit boundary i.e:
		[[z_0, z_1], # axial edit i=0
		[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
		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):
		"""
		Delegates methods called on this Core builder which do not match any method currently
		implemented in this builder to the contained Model instance.

		Args:
		attr: python method or attribute of the Model class
		"""
		return getattr(me.model, attr)


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

		@@ -73,7 +191,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].
		@@ -96,6 +214,10 @@ class MSRE(Core):
		assert(isinstance(dz, list))
		assert(len(dz)>0)
		assert(all(dz)>0.0)

		# call base class constructor
		super(MSRE, me).__init__(model, dz, sym_opt)

		if rodDomains:
		assert(utils.isArraySquare(rodDomains))
		assert(len(rodDomains)==len(coreMap))
		@@ -115,11 +237,6 @@ class MSRE(Core):


		thisMap = np.array(coreMap)

		me.dz = np.array(dz)

		me.model=model

		me.grBlockWidth = grBlockWidth

		# Use an equivalent diameter to get the same cross sectional area of the graphite block
		@@ -434,7 +551,6 @@ class MSRE(Core):
		chans = chansInDomain[domainID]
		model.addChannelsToDomain(domainID, chans)


		def plotGraphiteBlockLocations(fname):
		""" Make a plot of the core region.

		@@ -509,7 +625,6 @@ class MSRE(Core):
		obj.power = radPowers[i]
		i=i+1


		def getCoreEdgeBound(me):
		""" Returns the easternmost edge of the core assuming the origin is the center of the core map.

SubKit/build/Model.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -853,7 +853,7 @@ class Model:
		if chanVTK : me.editOptions['fluid_vtk'] = chanVTK
		if rodVTK : me.editOptions['rod_vtk'] = rodVTK
		if ctfHDF5 : me.editOptions['native_hdf5'] = ctfHDF5
		if veracsHDF5 : me.editOptions['hdf5'] = hdf5
		if veracsHDF5 : me.editOptions['hdf5'] = veracsHDF5
		if chanASCII : me.editOptions['chan_edits'] = chanASCII
		if rodEdits : me.editOptions['rod_edits'] = rodEdits
		if dnbEdits : me.editOptions['dnb_edits'] = dnbEdits

SubKit/build/Solid.py

+8 −11

Original line number	Diff line number	Diff line
		@@ -83,7 +83,6 @@ class FuelRod(SolidGeo):
		me.num_rings = num_rings
		me.percentTheoreticalDensity = percentTheoreticalDensity


		class _Material:
		""" Defines a solid material.

		@@ -158,5 +157,3 @@ class HeatedSolid(Solid):
		me.power = power
		me.powID = powID
		me.gapConductance = gapConductance

SubKit/build/SquareLatticeLWR_Nodal.py

+305 −12

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


		class NodalChannelGeom(object):
		""" Calculates geometry of a nodal channel in the model
		@@ -86,6 +88,7 @@ class NodalChannelGeom(object):
		""" Returns the size of nodal channel in y-direction"""
		return me.nodalChanSizes.getYsize(row, col)


		class NodalChannelBaseGeom(object):
		""" Calculates geometry of a nodal channel

		@@ -240,8 +243,10 @@ class NodalMesh(object):
		def getYLoc(me, row, col):
		return me.yLoc[row][col]


		class AssemblyChannelRelations(object):
		""" Determines relationship between nodal channels and the assemblies in the core
		This sometimes refered to as a chMapObj or CoreChanMap.

		Args:
		coreMap (list): 2D square list of integers. 0 means no assembly in the location. An integer
		@@ -249,7 +254,6 @@ class AssemblyChannelRelations(object):

		"""
		def __init__(me, coreMap):

		assert(utils.isArraySquare(coreMap))
		assert(utils.isSymmetric(coreMap))
		for row in coreMap:
		@@ -513,6 +517,7 @@ class NodalShroudData(object):
		else:
		return False


		class NodalChanDimensions(object):
		""" Calculates the x/y size of nodal channels in the model.

		@@ -609,11 +614,11 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		objects because they are not heated. Setting this to True will override the default of not modeling
		them. Note that whether they are modeled or not, their effect on gap size and wall friction are
		always captured in the model. Not modeling them will save some computational time.
		sym_opt (int): Enter either 1 for full core geometry or 4 for quarter symmetry geometry.
		"""
		def __init__(me, model, coreMap, assemMaps, pitch, dz, solidGeoms, startChanIndex=1, assemPitch=None,
		baffleGap=0.0, sectionID=1, symOpt=1, assemLosses=None, assemLossLevels=None, rodDomains=None,
		modelGuideTubes=False):

		modelGuideTubes=False, sym_opt=1):
		assert(isinstance(model, Model.Model))
		assert(len(coreMap)>0)
		assert(utils.isArraySquare(coreMap))
		@@ -623,6 +628,10 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		assert(len(dz)>0)
		assert(all(dz)>0.0)
		assert(startChanIndex>=1)
		assert(sym_opt==1)

		# call base class constructor
		super(SquareLatticeLWR_Nodal, me).__init__(model, dz, sym_opt)

		if isinstance(assemMaps, list):
		assemMapsDict = {1: assemMaps}
		@@ -717,8 +726,10 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		if isinstance(solidGeoDict[solidTypeID], Solid.FuelRod) or modelGuideTubes:
		me.rodsInAssem[assemID].append(solidID)
		[i, j]=chMapObj.getLoc(chIdx)
		model.addSolid(solidID, solidTypeID, Solid.HeatedSolid(mult=mult/4.0, power=power, powID=powID,
		x=nodalMesh.getXLoc(i, j), y=nodalMesh.getYLoc(i, j)))
		x, y = nodalMesh.getXLoc(i, j), nodalMesh.getYLoc(i, j)
		model.addSolid(solidID, solidTypeID,
		Solid.HeatedSolid(mult=mult/4.0, power=power, powID=powID,
		x=x, y=y))
		# Connect it ot the model
		model.addSolidOuterChan(solidID, chIdx, 1.0)

		@@ -793,11 +804,96 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		model.addRodsToDomain(domain, me.rodsInAssem[assemID], [domain]*len(me.rodsInAssem[assemID]))
		model.addChannelsToDomain(domain, chansInAssem[assemID])

		me.modelGuideTubes = modelGuideTubes
		model.addSection(sectionID, me.sec)
		me.model = model
		me.chMapObj = chMapObj
		me.nodalGeom = nodalGeom
		me.dz = dz

		def setEditOptions(me, chanVTK=None, rodVTK=None, ctfHDF5=None, veracsHDF5=None, chanASCII=None,
		rodEdits=None, dnbEdits=None):
		"""
		Delegates setting edits to the Model instance contained in this builder.

		Args:
		chanVTK (bool): Set to True to write the channel VTK file (default True).
		rodVTK (bool): Set to True to write the rod VTK file (default False).
		ctfHDF5 (bool): Set to True to write the native CTF HDF5 file (default True).
		veracsHDF5 (bool): Set to True to write the VERA-CS HDF5 file
		(requires that core map be specified and will only work for square rod lattices) (default False).
		chanASCII (bool): Set to True to write the ASCII channel output file.
		rodEdits (bool): Set to True to write rod edit information to the standard CTF .out file
		(Note the rod information in the .out file is limited in scope and it is better to
		get rod data from the native HDF5 file).
		dnbEdits (bool): Set to True to produce the CTF .dnb file, which contains information about DNB (default False).
		"""
		if veracsHDF5 == True:
		rodMap, chanMap, assemMap, symOpt = me.getChannelAndRodMaps()
		me.model.setCoreMaps(rodMap, chanMap, assemMap, symOpt)
		if me.modelGuideTubes:
		raise RuntimeError("Nodal VERA CS edits are not compatible with modelGuideTubes=True. " + \
		" Set modelGuideTubes=False")
		me.model.setEditOptions(chanVTK, rodVTK, ctfHDF5, veracsHDF5, chanASCII, rodEdits, dnbEdits)

		def getChannelAndRodMaps(me):
		"""
		Returns the rodMap and chanMap required by group 17 for veracsHDF5 file
		writing.

		Example use:
		>>> my_lwr_nodal = SquareLatticeLWR_Nodal(...)
		>>> rodMap, chanMap, assemMap, symOpt = my_lwr_nodal.getChannelAndRodMaps()
		>>> my_lwr_nodal.model.setCoreMaps(rodMap, chanMap, assemMap, symOpt)
		>>> my_lwr_nodal.model.setEditOptions(veracsHDF5=True)

		Returns:
		rodMap, chanMap, assemMap, symOpt
		rodMap is a 2 dimensional list
		chanMap is a 2 dimensional list
		assemMap is a 2 dimensional list
		symOpt is an int: 1=full sym, 4=qtr sym
		"""
		if me.modelGuideTubes:
		raise RuntimeError("Model guide tubes is currently not supported " + \
		"by the nodal veracsHDF5 writer. Set modelGuideTubes=False.")
		else:
		rodMap = np.zeros((len(me.chMapObj.chAssemIdx),len(me.chMapObj.chAssemIdx[0])), dtype=int)

		assert me.core_sym == 1
		for obj in me.model.solidObjects.values():
		# this only works in full sym mode for now
		assert(obj.symtype==1)

		# create the rodMap
		for asm_i in range(me.chMapObj.getAssemDimLen()):
		for asm_j in range(me.chMapObj.getAssemDimLen()):
		assemID = me.chMapObj.getAssemIndexFromLoc(asm_i, asm_j)
		if assemID:
		rods = me.rodsInAssem[assemID]
		if me.modelGuideTubes: assert len(rods) == 8
		if not me.modelGuideTubes: assert len(rods) == 4

		# fill rodMap
		if me.modelGuideTubes:
		# top nodes in asm
		rodMap[asm_i * 2 + 0][asm_j * 4 + 0] = int(assemID)
		rodMap[asm_i * 2 + 0][asm_j * 4 + 1] = int(assemID)
		rodMap[asm_i * 2 + 0][asm_j * 4 + 2] = int(assemID)
		rodMap[asm_i * 2 + 0][asm_j * 4 + 3] = int(assemID)
		# bottom nodes in asm
		rodMap[asm_i * 2 + 1][asm_j * 4 + 0] = int(assemID)
		rodMap[asm_i * 2 + 1][asm_j * 4 + 1] = int(assemID)
		rodMap[asm_i * 2 + 1][asm_j * 4 + 2] = int(assemID)
		rodMap[asm_i * 2 + 1][asm_j * 4 + 3] = int(assemID)
		else:
		# top nodes in asm
		rodMap[asm_i * 2 + 0][asm_j * 2 + 0] = int(assemID)
		rodMap[asm_i * 2 + 0][asm_j * 2 + 1] = int(assemID)
		# bottom nodes in asm
		rodMap[asm_i * 2 + 1][asm_j * 2 + 0] = int(assemID)
		rodMap[asm_i * 2 + 1][asm_j * 2 + 1] = int(assemID)
		# convert rodMap to list of lists to be consistent with SubKit
		rodMap = list([list(rodRow) for rodRow in rodMap])
		return rodMap, me.chMapObj.chAssemIdx, me.chMapObj.assemIdx, me.core_sym

		def setCoreInletBC(me, massflux, temperature, flowRamp=None):
		""" Sets the mass flux and temperature uniformly at the inlet of the core.
		@@ -831,15 +927,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
		@@ -893,3 +989,200 @@ 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):
		"""
		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)

		An example node layout for a 4 assembly model is:
		-------------
		\| 1 2 \| 1 2 \|
		\| 3 4 \| 3 4 \|
		-------------
		\| 1 2 \| 1 2 \|
		\| 3 4 \| 3 4 \|
		-------------
		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))

		# Renormalize the radial powers
		radP = []
		mults = []
		for obj in me.model.solidObjects.values():
		radP.append(obj.power)
		mults.append(obj.mult)
		radP = np.array(radP)/np.average(np.array(radP), weights=mults)
		for i, obj in enumerate(me.model.solidObjects.values()):
		obj.power = radP[i]

		def _setCorePowerShapeAsm3D(me, powerDist, wgts=None, coreStart=0.0):
		"""
		Set a 3D power distribution for the core. Each assembly can have its own
		axial power shape.

		Args:
		powerDist (numpy.ndarray): 3 dim array. Contains power distribution with indexing
		(asm_i, asm_j, 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]
		"""
		assert len(powerDist.shape) == 3
		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_centers.size == powerDist.shape[2]

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

		radPowers = []
		weights = []
		for obj in me.model.solidObjects.values():
		assert(obj.symtype==1)
		for i in range(me.chMapObj.getAssemDimLen()):
		for j in range(me.chMapObj.getAssemDimLen()):
		assemID = me.chMapObj.getAssemIndexFromLoc(i, j)
		if assemID:
		powID = "pwr_" + str(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
		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_axpwr

		def _setCorePowerShapeNodal3D(me, powerDist, wgts=None, coreStart=0.0):
		"""
		Set a 3D power distribution for the core. Each node within each assembly
		can have its own axial power shape.

		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:
		intra-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
		so the power tables are correctly defined in CTF. [m]
		"""
		assert len(powerDist.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]

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

		radPowers = []
		weights = []
		for obj in me.model.solidObjects.values():
		assert(obj.symtype==1)
		for asm_i in range(me.chMapObj.getAssemDimLen()):
		for asm_j in range(me.chMapObj.getAssemDimLen()):
		assemID = me.chMapObj.getAssemIndexFromLoc(asm_i, asm_j)
		if assemID:
		rods = me.rodsInAssem[assemID]
		if me.modelGuideTubes: assert len(rods) == 8
		if not me.modelGuideTubes: assert len(rods) == 4

		# determine quadrent of rod in assem
		rod_xy = {}
		for rodID in rods:
		obj_rod = me.model.solidObjects[rodID]
		rod_xy[rodID] = (obj_rod.x, obj_rod.y)
		unique_xy = np.unique(np.asarray(rod_xy.values()), axis=0)
		# 4 pairs of unique (x,y) coords must be present
		assert unique_xy.size == 4 * 2

		centroid_xy = np.average(unique_xy, axis=0)
		for rodID in rods:
		obj_rod = me.model.solidObjects[rodID]
		if obj_rod.x > centroid_xy[0] and obj_rod.y > centroid_xy[1]:
		# north east
		nodeID = 1
		elif obj_rod.x < centroid_xy[0] and obj_rod.y > centroid_xy[1]:
		# north west
		nodeID = 0
		elif obj_rod.x < centroid_xy[0] and obj_rod.y < centroid_xy[1]:
		# south west
		nodeID = 2
		else:
		# south east
		nodeID = 3

		# set rod power
		powID = "pwr_asm_" + str(assemID) + "_node_" + str(nodeID)
		node_axial_power_profile = powerDist[asm_i, asm_j, nodeID, :]

		# 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_axpwr

SubKit/build/genDeck.py

+9 −2

File changed.

Preview size limit exceeded, changes collapsed.