build chanMap and rodMap for veracsHDF5 output (03b8380d) · Commits · Phathanapirom, Birdy / SubKit

SubKit/build/SquareLatticeLWR_Nodal.py

+94 −7

Original line number	Diff line number	Diff line
		@@ -241,8 +241,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
		@@ -250,7 +252,6 @@ class AssemblyChannelRelations(object):

		"""
		def __init__(me, coreMap):

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


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

		@@ -613,8 +615,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		"""
		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, setCoreMaps=False):
		assert(isinstance(model, Model.Model))
		assert(len(coreMap)>0)
		assert(utils.isArraySquare(coreMap))
		@@ -796,11 +797,87 @@ 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
		if setCoreMaps:
		rodMap, chanMap, assemMap, symOpt = me.getChannelAndRodMaps()
		me.model.setCoreMaps(rodMap, chanMap, assemMap, symOpt)

		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:
		where rodMap is a 2 dimensional list
		where chanMap is a 2 dimensional list
		where assemMap is a 2 dimensional list
		symOpt is an int: 1=full sym, 4=qtr sym
		"""
		if not me.modelGuideTubes:
		rodMap = np.zeros((len(me.chMapObj.chAssemIdx),len(me.chMapObj.chAssemIdx[0])), dtype=int)
		else:
		rodMap = np.zeros((len(me.chMapObj.chAssemIdx),2*len(me.chMapObj.chAssemIdx[0])), dtype=int)
		symOpt = 1
		for obj in me.model.solidObjects.values():
		# this only works in full sym mode for now
		assert(obj.symtype==1)
		symOpt = obj.symtype

		# 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

		tmp_rod_types = []
		tmp_rod_x = []
		tmp_rod_y = []
		for rodID in rods:
		obj_rod = me.model.solidObjects[rodID]
		obj_rod_type_id = me.model.solidTypeIDs[rodID]
		tmp_rod_types.append(obj_rod_type_id)
		tmp_rod_x.append(obj_rod.x)
		# in nodal model: xy loc of guide tube == xy locs of heated rods
		tmp_rod_y.append(obj_rod.y)

		# 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(0)
		rodMap[asm_i * 2 + 0][asm_j * 4 + 2] = int(assemID)
		rodMap[asm_i * 2 + 0][asm_j * 4 + 3] = int(0)
		# 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(0)
		rodMap[asm_i * 2 + 1][asm_j * 4 + 2] = int(assemID)
		rodMap[asm_i * 2 + 1][asm_j * 4 + 3] = int(0)
		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, symOpt

		def setCoreInletBC(me, massflux, temperature, flowRamp=None):
		""" Sets the mass flux and temperature uniformly at the inlet of the core.
		@@ -958,11 +1035,11 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):

		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 modelGuideTubes=True)
		# obj_rod_type_id == 2 is guide tube (only present if me.modelGuideTubes=True)
		# obj_rod is type SubKit.build.Solid.HeatedSolid
		obj_rod.powID = powID
		# radial power factors are implicit in the 3d power distribution
		obj_rod.power = 1.0

		@@ -1011,6 +1088,8 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		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 = {}
		@@ -1018,11 +1097,14 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		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)
		assert unique_xy.size == 4 * 2

		# guide tubes are offset from heated rods in nodal mesh
		if me.modelGuideTubes: assert unique_xy.size == 4 * 2
		if not me.modelGuideTubes: assert unique_xy.size == 4

		centroid_xy = np.average(unique_xy, axis=0)
		for rodID in rods:
		obj_rod = me.model.solidObjects[rodID]
		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 = 1
		@@ -1041,4 +1123,9 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		me.model.addAxialPowerShape(powID, me.z_power_grid,
		powerDist[asm_i, asm_j, nodeID, :])
		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 = 1.0

SubKit/build/genDeck.py

+5 −1

Original line number	Diff line number	Diff line
		@@ -97,7 +97,11 @@ def writeDeck(model, filename):
		fluidprops=11
		else:
		assert(False)
		group1Data.append(" 1 2 0 3{:5d}{:15.5e}{:6d} 1 0{:5d}{:5d} 0 7 0\n".format(model.solver, mdotInit, notrans, fluidprops, mflx))
		if model.rodMap:
		hasmaps = 1
		else:
		hasmaps = 0
		group1Data.append(" 1 2 0 3{:5d}{:15.5e}{:6d} 1{:5d}{:5d}{:5d} 0 7 0\n".format(model.solver, mdotInit, notrans, hasmaps, fluidprops, mflx))
		group1Data.append("*Card 1.2\n")
		if mflx==1:
		group1Data.append("** GTOT AFLUX DHFRAC MASSFLUX\n")

tests/nodal_p7/buildDeck.py

+1 −0

Original line number	Diff line number	Diff line
		@@ -74,6 +74,7 @@ def main():
		[ 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
		[ 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]]


		# Fuel rod is type 1, but mark it with an "f" to make it more visible
		# Guide tube is type 2
		f=1

tests/nodal_p7_pow3d/buildDeck.py

+18 −2

Original line number	Diff line number	Diff line
		@@ -189,6 +189,22 @@ def main(power_3d_flag=1):
		# [ 0, 0, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 0, 0],
		# [ 0, 0, 0, 0, 3, 3, 3, 4, 4, 4, 4, 0, 0, 0, 0]]

		# crate assembly map for vera HDF5 edits
		asmMap = np.zeros((15,15), dtype=int)
		n_asm = 0
		for i in range(asmMap.shape[0]):
		for j in range(asmMap.shape[1]):
		if coreMap[i][j] >= 1:
		n_asm += 1
		asmMap[i, j] = n_asm
		asmMap = [list(asmRow) for asmRow in asmMap]

		# crate rod map for vera HDF5 edits
		pass

		# crate channel map for vera HDF5 edits
		pass

		print np.unique(np.array(rodDomains), return_counts=True)

		model = cdb.Model()
		@@ -232,7 +248,7 @@ def main(power_3d_flag=1):

		builder = SquareLatticeLWR_Nodal(model, coreMap, pwr17x17, pitch, dz, solidGeoms, assemPitch=assemPitch,
		baffleGap=baffleGap, assemLosses=formLoss, assemLossLevels=gridLevels, rodDomains=rodDomains,
		modelGuideTubes=True)
		modelGuideTubes=True, setCoreMaps=True)

		# Run to 1 second to let the solution reach steady state
		model.addTransientGroup(tEnd=1.0)
		@@ -254,7 +270,7 @@ def main(power_3d_flag=1):
		model.setAveragePower(qprime=linearHeatRate, dhfrac=directHeat)

		# Turn on rod VTK and DNB edits. They are both off by default.
		model.setEditOptions(dnbEdits=True,veracsHDF5=True)
		model.setEditOptions(dnbEdits=True,veracsHDF5=True,ctfHDF5=True)

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