Merge remote-tracking branch 'willfork/set_3dpower_dist' into calc_nodal_powers (9da3f37f) · Commits · Phathanapirom, Birdy / SubKit

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

+148 −2

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


		class NodalChannelGeom(object):
		""" Calculates geometry of a nodal channel in the model

		@@ -717,8 +718,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)

		@@ -893,3 +896,146 @@ 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. Each assembly can have it's 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_power_grid.size == powerDist.shape[2]

		# normalize 3d powers
		if wgts is None:
		wgts = np.ones(powerDist.shape)
		assert wgts.shape == powerDist.shape
		powerDist *= wgts
		powerDist /= np.average(powerDist, weights=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)
		me.model.addAxialPowerShape(powID, me.z_power_grid, powerDist[i, j, :])
		rods = me.rodsInAssem[assemID]

		for rodID in rods:
		obj_rod = me.model.solidObjects[rodID]
		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 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

		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.

		Args:
		powerDist (numpy.ndarray): 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]
		"""
		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_power_grid.size == powerDist.shape[2]

		# normalize 3d powers
		if wgts is None:
		wgts = np.ones(powerDist.shape)
		assert wgts.shape == powerDist.shape
		powerDist *= wgts
		powerDist /= np.average(powerDist, weights=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]

		# 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)
		assert unique_xy.size == 4 * 2
		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 = 0
		elif obj_rod.x < centroid_xy[0] and obj_rod.y > centroid_xy[1]:
		# north west
		nodeID = 1
		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)
		me.model.addAxialPowerShape(powID, me.z_power_grid,
		powerDist[asm_i, asm_j, nodeID, :])
		obj_rod.powID = powID
		obj_rod.power = 1.0

SubKit/process/Hdf5Tools.py

+1 −2

Original line number	Diff line number	Diff line
		@@ -273,7 +273,6 @@ class Hdf5Interface(object):
		me._pinIndexValid(pin)
		return me.pinAxialNodes[0:me.pinNumAxial[pin-1], pin-1]


		def getChanArea(me, ch):
		""" Returns the nominal channel flow area [m**2]

		@@ -282,7 +281,7 @@ class Hdf5Interface(object):

		"""
		me._chanIndexValid(ch)
		return me.h5['CORE/chan_area'][ch-1].value
		return me.h5['CORE/chan_area'][ch-1]

		def getPinThetas(me, pin, level):
		""" Returns the azimuthal fractions for the passed pin and level

tests/nodal_p7_pow3d/buildDeck.py

0 → 100644

+283 −0

Original line number	Diff line number	Diff line
		#!/usr/bin/env python
		# Author: Bob Salko
		# Date: 1/23/19
		# Description:
		# Makes a model of CASL problem 7 using 4 subchannels per assembly
		#
		import numpy as np
		import math
		import copy
		from SubKit.build.SquareLatticeLWR_Nodal import SquareLatticeLWR_Nodal
		import SubKit.utils.UnitConversions as units
		import SubKit.build as cdb
		import sys
		sys.path.insert(0, '../../')


		def main(power_3d_flag=True):

		dFuelOuter = 0.475e-2*2 # m
		dFuelInner = 0.418e-2*2 # m
		dFuelPellet = 0.4096e-2*2 # m
		dGuideTubeOuter = 0.602e-2*2 # m
		dGuideTubeInner = 0.561e-2*2 # m
		assemPitch = 21.5e-2 # m
		pitch = 1.26e-2 # m
		baffleGap = 0.19e-2 # m
		inletMassFlux = 0.3483429E+04 # kg/m**2/s
		inletTemp = 291.85 # C
		outletPressure = 155.132039 # bar
		linearHeatRate = 16.71977 # kW/m
		directHeat = 0.026
		gridLocations = list(
		np.array([13.884, 75.2, 127.4, 179.6, 231.8, 284.0, 336.2])*units.t_cm_m)
		gridHeights = list(
		np.array([3.866, 3.810, 3.810, 3.810, 3.810, 3.810, 3.810])*units.t_cm_m)
		axial_edit_bounds = list(np.array([11.951, 15.817, 24.028, 32.239, 40.45, 48.662, 56.873, 65.084, 73.295, 77.105, 85.17,
		93.235, 101.3, 109.365, 117.43, 125.495, 129.305, 137.37, 145.435, 153.5, 161.565, 169.63, 177.695,
		181.505, 189.57, 197.635, 205.7, 213.765, 221.83, 229.895, 233.705, 241.77, 249.835, 257.9, 265.965, 274.03,
		282.095, 285.905, 293.97, 302.035, 310.1, 318.165, 326.23, 334.295, 338.105, 346.0262, 353.9474, 361.8686,
		369.7898, 377.711])*units.t_cm_m)
		formLoss = 0.9070

		# Setup axial mesh
		gridLocations = list(np.array(gridLocations)-axial_edit_bounds[0])
		axial_edit_bounds = list(np.array(axial_edit_bounds)-axial_edit_bounds[0])
		dz = []
		for j in range(1, len(axial_edit_bounds)):
		dz.append(axial_edit_bounds[j]-axial_edit_bounds[j-1])
		numLevels = len(dz)
		gridBottoms = list(np.array(gridLocations)-np.array(gridHeights)/2.0)
		gridLevels = []
		for grid in gridBottoms:
		found = False
		for level, z in enumerate(axial_edit_bounds):
		if np.isclose(grid, z):
		gridLevels.append(level+1)
		found = True
		assert(found)
		assert(len(gridLevels) == len(gridLocations))

		# The map of the core. Each 1 represents an assembly
		coreMap = [[0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
		[0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
		[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
		[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
		[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
		[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
		[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
		[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
		[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
		[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
		[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
		[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
		[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
		[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
		pwr17x17 = [
		[f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
		[f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
		[f, f, f, f, f, 2, f, f, 2, f, f, 2, f, f, f, f, f],
		[f, f, f, 2, f, f, f, f, f, f, f, f, f, 2, f, f, f],
		[f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
		[f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f],
		[f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
		[f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
		[f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f],
		[f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
		[f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
		[f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f],
		[f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
		[f, f, f, 2, f, f, f, f, f, f, f, f, f, 2, f, f, f],
		[f, f, f, f, f, 2, f, f, 2, f, f, 2, f, f, f, f, f],
		[f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f],
		[f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f]]

		# The multiple gap connects are not working with CTF parallel for some reason
		rodDomains = [[0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0],
		[0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0],
		[0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
		[0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
		[1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
		[1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
		[1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
		[3, 3, 3, 3, 3, 3, 3, 4, 2, 2, 2, 2, 2, 2, 2],
		[3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4],
		[3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4],
		[3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4],
		[0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0],
		[0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0],
		[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]]
		# rodDomains = [[ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0],
		# [ 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0],
		# [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
		# [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
		# [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
		# [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
		# [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
		# [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
		# [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
		# [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
		# [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2],
		# [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
		# [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
		# [ 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0],
		# [ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0]]
		# rodDomains = [[ 0, 0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0],
		# [ 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 0, 0],
		# [ 0, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 0],
		# [ 0, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 0],
		# [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3],
		# [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3],
		# [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3],
		# [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 6, 6, 6, 6, 6, 6],
		# [ 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6],
		# [ 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6],
		# [ 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6],
		# [ 0, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 0],
		# [ 0, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 0],
		# [ 0, 0, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 0, 0],
		# [ 0, 0, 0, 0, 4, 4, 5, 5, 5, 6, 6, 0, 0, 0, 0]]
		# rodDomains = [[ 0, 0, 0, 0, 1, 2, 2, 2, 2, 3, 3, 0, 0, 0, 0],
		# [ 0, 0, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 0, 0],
		# [ 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 0],
		# [ 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 0],
		# [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3],
		# [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3],
		# [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3],
		# [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 6, 6, 6, 6, 6],
		# [ 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6],
		# [ 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6],
		# [ 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6],
		# [ 0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0],
		# [ 0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0],
		# [ 0, 0, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 0, 0],
		# [ 0, 0, 0, 0, 4, 5, 5, 5, 5, 5, 6, 0, 0, 0, 0]]
		rodDomains = [[0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 3, 0, 0, 0, 0],
		[0, 0, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 0, 0],
		[0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 0],
		[0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 0],
		[1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3],
		[1, 1, 1, 1, 7, 7, 7, 2, 2, 2, 3, 3, 3, 3, 3],
		[1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 3],
		[1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6],
		[4, 4, 4, 4, 4, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6],
		[4, 4, 4, 4, 4, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6],
		[4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6],
		[0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0],
		[0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0],
		[0, 0, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 0, 0],
		[0, 0, 0, 0, 4, 5, 5, 5, 5, 5, 6, 0, 0, 0, 0]]
		# rodDomains = [[ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0],
		# [ 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0],
		# [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
		# [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0],
		# [ 1, 1, 1, 1, 1, 1, 5, 5, 5, 2, 2, 2, 2, 2, 2],
		# [ 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 2, 2, 2, 2, 2],
		# [ 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 2, 2, 2, 2],
		# [ 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 2, 2, 2, 2],
		# [ 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4],
		# [ 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4],
		# [ 3, 3, 3, 3, 3, 3, 5, 5, 5, 4, 4, 4, 4, 4, 4],
		# [ 0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0],
		# [ 0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0],
		# [ 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]]

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

		model = cdb.Model()

		model.setTitle("4-Loop nodal mesh loss of flow")

		solidGeoms = {1: cdb.FuelRod(d_outer=dFuelOuter, d_inner=dFuelInner, d_pellet=dFuelPellet), # Fuel rod
		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],
		[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,
		0.97, 1.17, 0.93, 1.12, 0.74, 0.00, 0.00],
		[0.00, 0.75, 1.12, 0.93, 1.12, 0.94, 1.19, 0.93,
		1.19, 0.94, 1.12, 0.93, 1.12, 0.75, 0.00],
		[0.75, 1.10, 0.93, 1.12, 0.95, 1.19, 0.95, 1.22,
		0.95, 1.19, 0.95, 1.12, 0.93, 1.10, 0.75],
		[0.76, 0.93, 1.17, 0.94, 1.19, 0.98, 1.19, 0.94,
		1.19, 0.98, 1.19, 0.94, 1.17, 0.93, 0.76],
		[0.77, 1.10, 0.97, 1.19, 0.95, 1.19, 0.94, 1.12,
		0.94, 1.19, 0.95, 1.19, 0.97, 1.10, 0.77],
		[0.78, 0.91, 1.18, 0.93, 1.22, 0.94, 1.12, 0.90,
		1.12, 0.94, 1.22, 0.93, 1.18, 0.91, 0.78],
		[0.77, 1.10, 0.97, 1.19, 0.95, 1.19, 0.94, 1.12,
		0.94, 1.19, 0.95, 1.19, 0.97, 1.10, 0.77],
		[0.76, 0.93, 1.17, 0.94, 1.19, 0.98, 1.19, 0.94,
		1.19, 0.98, 1.19, 0.94, 1.17, 0.93, 0.76],
		[0.75, 1.10, 0.93, 1.12, 0.95, 1.19, 0.95, 1.22,
		0.95, 1.19, 0.95, 1.12, 0.93, 1.10, 0.75],
		[0.00, 0.75, 1.12, 0.93, 1.12, 0.94, 1.19, 0.93,
		1.19, 0.94, 1.12, 0.93, 1.12, 0.75, 0.00],
		[0.00, 0.00, 0.74, 1.12, 0.93, 1.17, 0.97, 1.18,
		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]]

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

		# Run to 1 second to let the solution reach steady state
		model.addTransientGroup(tEnd=1.0)
		# Run a 10 second transient
		model.addTransientGroup(tEnd=11.0, editInterval=1.0)

		# Ramp the power
		time = [1.0, 4.0]
		power = [1.0, 1.0]
		model.addTotalPowerRamp(time, power)

		# Set the core boundary conditions
		time = [1.0, 1.1]
		flow = [1.00, 1.00]
		builder.setCoreInletBC(inletMassFlux, inletTemp, flowRamp=[time, flow])
		builder.setCoreOutletBC(outletPressure, inletTemp)

		# Set the nominal power and direct heating to coolant
		model.setAveragePower(qprime=linearHeatRate, dhfrac=directHeat)

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

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

		def axPow(z):
		return np.cos((z-coreHeight/2)*np.pi/2/coreHeight)

		if power_3d_flag:
		# set 3d power distribution
		axialMesh = builder.z_centers
		radPow_arr = np.array(radPow)
		radPow_arr_full = np.ones(
		(radPow_arr.shape[0], radPow_arr.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
		builder.setCorePowerShape3D(powerDist3d)
		model.generateModel('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)

tests/unitTests/test_Hdf5Tools.py

+8 −0

Original line number	Diff line number	Diff line
		@@ -485,6 +485,14 @@ class test_Hdf5Tools(unittest.TestCase):
		Tave = np.average(Temp, weights=mdotTot)
		me.assertAlmostEqual(Tave, me.obj.getAveTemp(section=2, level=2, state=2))

		def test_getChanArea(me):
		chArea = [8.536799999999999E-5, 8.536799999999999E-5, 8.536799999999999E-5,
		8.536799999999999E-5, 8.536799999999999E-5, 8.536799999999999E-5,
		8.536799999999999E-5, 1.7073999999999998E-4, 1.7073999999999998E-4,
		1.7073999999999998E-4]
		for i, ch in enumerate(range(10)):
		me.assertAlmostEqual(me.obj.getChanArea(ch + 1), chArea[i])

		def test_getChanVoid(me):
		axial, void = me.obj.getChanMeshAndData('chan_void_vap', chan=2, state=2)
		state2_chan2_level4 = 1.0007418278303053E-6