Several improvements for the nodal model builder (f5a55f53) · Commits · Gurecky, William / SubKit

SubKit/build/Model.py

+58 −6

Original line number	Diff line number	Diff line
		@@ -101,9 +101,19 @@ class Model:

		me.modelParallel = False

		# Stopping criteria for steady state
		# Define a steady state model and write time group data and
		# stopping criteria appropriate for steady state
		me.runSteadyState = True

		# Steady-state convergence criteria are not normally printed
		# for transient models and are normally ignored by CTF for
		# transient models. However, for certain applications
		# (e.g. code coupling) the convergence criteria may still be
		# used even if the input file specifies the model as a
		# transient. This flag forces printing of the convergence
		# criteria even for transient models.
		me.transPrintConv = False

		# Is the model nodal
		me._is_nodal = False

		@@ -546,6 +556,39 @@ class Model:
		me.qprime = qprime
		me.dhfrac = dhfrac

		def setTotalCorePower(me, power, dhfrac=0.0):
		""" Takes the total core power and sets the average linear heat rate per rod in the model.

		Args:
		power (float): Total core power [kW]
		dhfrac (float): The direct heating fraction. Must be between 0 and 1. This percentage of
		heat generated in the rods will be directly deposited into the fluid and bypass going through
		the rod surface.
		"""
		nRods = me.getTotalFuelRods()
		if nRods == 0:
		raise RuntimeError('setTotalCorePower must be called after the solid geometries are defined.')
		coreHeight = me._getModelHeight()
		qprime = power / (coreHeight * nRods)
		me.setAveragePower(qprime, dhfrac)

		def getTotalPins(me):
		""" Returns the total number of heated pins (including pin multipliers) in the model """
		nPins = 0.
		for (ipin,pin) in me.solidObjects.items():
		if isinstance(pin, Solid.HeatedSolid):
		nPins += pin.mult
		return nPins

		def getTotalFuelRods(me):
		""" Returns the total number of fuel rods (including pin multipliers) in the model """
		nRods = 0.
		for (ipin,pin) in me.solidObjects.items():
		objType = me.solidTypeIDs[ipin]
		if isinstance(me.solidTypes[objType], Solid.FuelRod):
		nRods += pin.mult
		return nRods

		def addSection(me, secID, section):
		""" Use this to add a section to the model.

		@@ -1199,7 +1242,7 @@ class Model:
		assert gapID not in me.gapGeoTables
		me.gapGeoTables[gapID] = tableID

		def addTransientGroup(me, tEnd, editInterval=None, dtMin=None, dtMax=None):
		def addTransientGroup(me, tEnd, editInterval=None, dtMin=None, dtMax=None, rtwfp=1.0):
		"""
		Adds a transient group to the simulation

		@@ -1216,6 +1259,10 @@ class Model:
		editInterval (float): The time between code edits (defaults to tEnd) [s]
		dtMin (float): Minimum allowable timestep size for this group [s]
		dtMax (float): Maximum allowable timestep size for this group [s]
		rtwfp (float): Ratio of timestep size between solid and fluid solutions. Note: under
		normal circumstances, CTF will reset rtwfp to 1 for transients. Setting
		a value other than 1 will only impact the calculation under specific
		circumstances (e.g. system-coupled simulations)

		"""
		me.runSteadyState = False
		@@ -1234,7 +1281,7 @@ class Model:
		dtMax = me.dtMax_DEFAULT
		me.dtMin.append(dtMin)
		me.dtMax.append(dtMax)
		me.rtwfp = 1.0
		me.rtwfp = rtwfp

		def setSolver(me, solver):
		""" Sets the solver to use for solving the pressure matrix.
		@@ -1264,7 +1311,7 @@ class Model:
		assert maxOuterIterations > 0
		me.numericalControls['maxOuterIterations'] = maxOuterIterations

		def setSteadyState(me, energyBalance=None, massBalance=None, maxIterations=None, rtwfp=None,
		def setSteadyState(me, transPrintConv=False, energyBalance=None, massBalance=None, maxIterations=None, rtwfp=None,
		p_l2=None, p_linf=None, Tcool_l2=None, Tcool_linf=None,
		Tsolid_l2=None, Tsolid_linf=None, vl_l2=None, vl_linf=None, vv_l2=None, vv_linf=None,
		vd_l2=None, vd_linf=None, void_l2_abs=None, void_linf_abs=None, p_l2_abs=None, p_linf_abs=None, Tcool_l2_abs=None,
		@@ -1273,7 +1320,8 @@ class Model:
		dtmin=None, dtmax=None):
		""" Call to specify the model as steady state.

		Throws exception if transient groups specified already.
		Throws exception if transient groups specified already. This exception can be bypassed using
		transPrintConv, which allows card group 19 criteria to be printed even for a transient deck.
		All stopping criteria have defaults that will be used if not explicitly set in the argument list.
		The stopping criteria on solution variables (e.g., pressure, velocity, temperature) are all unitless.
		The criteria ensure that the solution variables in CTF are changing less the specified amount between
		@@ -1281,6 +1329,8 @@ class Model:
		a single value. The CTF Theory Manual can be reviewed for details on how the metrics are formulated.

		Args:
		transPrintConv (bool): Flag to allow printing of card group 19 regardless of whether this is a
		steady state or transient deck
		energyBalance (float): Tolerance on energy balance (energy in/energy out) [unitless]
		massBalance (float): Tolerance on mmass balance (mass in/mass out) [unitless]
		maxIterations (int): Maximum iterations CTF should take before stopping the simulation [unitless]
		@@ -1316,7 +1366,9 @@ class Model:
		dtmax (float): Maximum allowable timestep size (steady-state is modeled as a pseudo-transient)

		"""
		if not me.runSteadyState:
		me.transPrintConv = transPrintConv

		if not me.runSteadyState and not transPrintConv:
		raise RuntimeError(
		'Already called addTransientGroup, so cannot convert to steady state run')

SubKit/build/SquareLatticeLWR_Nodal.py

+47 −26

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


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

		@@ -662,13 +661,22 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		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.
		gapWidthMult (float): Global multiplier that will be applied to all gap width values in the model.
		gapWidthMultNarrow (float): Global multiplier that will be applied to all gaps with a nominal width
		less than or equal to gapWidthMultThreshold. Useful for applying a multiplier specifically to
		intra-assembly gaps (which are typically < 0.01 m)
		gapWidthMultWide (float): Global multiplier that will be applied to all gaps with a nominal width
		greater than gapWidthMultThreshold. Useful for applying a multiplier specifically to inter-assembly
		gaps (which are typically > 0.01 m)
		gapWidthMultThreshold (float): Threshold for applying narrow or wide gap width multipliers [m]. Each gap with
		a nominal width less than gapWidthMultThreshold will be multiplied by gapWidthMultNarrow, otherwise
		it will be multiplied by gapWidthMultWide. Defaults to 0.01 m, which for typical 17x17 PWR fuel should
		successfully distinguish between all intra- and inter-assembly gaps.
		gapLengthMult (float): Global multiplier that will be applied to all gap lengths in the model.
		"""
		def __init__(me, model, coreMap, assemMaps, pitch, dz, solidGeoms, startChanIndex=1, assemPitch=None,
		baffleGap=0.0, sectionID=1, symOpt=1, assemLosses=None, assemLossLevels=None,
		assemGridBlockageRatios=None, assemGridEnhancementCoeffs=None, rodDomains=None, modelGuideTubes=False, sym_opt=1,
		gapWidthMult=1.0, gapLengthMult=1.0):
		gapWidthMultNarrow=1.0, gapWidthMultWide=1.0, gapLengthMult=1.0, gapWidthMultThreshold=0.01):
		assert(isinstance(model, Model.Model))
		assert(len(coreMap) > 0)
		assert(utils.isArraySquare(coreMap))
		@@ -740,7 +748,11 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		if southGaps:
		gaps = gaps + southGaps
		for gap in gaps:
		model.setGap(gap.ii, gap.jj, gap.width*gapWidthMult,
		if gap.width > gapWidthMultThreshold:
		width = gap.width * gapWidthMultWide
		else:
		width = gap.width * gapWidthMultNarrow
		model.setGap(gap.ii, gap.jj, width,
		gap.length*gapLengthMult, mult=gap.mult)

		# Add the solid geometry objects
		@@ -776,7 +788,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		power = 1.0
		elif modelGuideTubes:
		powID = 0
		power = 0.0
		power = 1.0
		if isinstance(solidGeoDict[solidTypeID], Solid.FuelRod) or modelGuideTubes:
		me.rodsInAssem[assemID].append(solidID)
		[i, j] = chMapObj.getLoc(chIdx)
		@@ -1128,6 +1140,31 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		me.model.setInitialConditions(
		me.mdotTotal, me.tempInitial, me.outletPressure)

		def normalizeRadialPowers(me):
		# Renormalize the radial powers over all fuel rods. Set the
		# radial power factor for unheated pins to 1.0 (these pins
		# will be assigned a zero axial power shape, so their power
		# will still come out zero).
		radP = []
		powMults = []
		for i, obj in enumerate(me.model.solidObjects.values()):
		solidID = me.model.solidObjects.keys()[i]
		radP.append(obj.power)
		objType = me.model.solidTypeIDs[solidID]
		# powMults contains the number of fuel rod pins in each
		# solid (does not include unheated pins like obj.mult
		# does)
		if isinstance(me.model.solidTypes[objType], Solid.FuelRod):
		powMults.append(obj.mult)
		else:
		powMults.append(0.)
		radP = np.array(radP) / np.average(np.array(radP), weights=powMults)
		for i, obj in enumerate(me.model.solidObjects.values()):
		if powMults[i] > 0:
		obj.power = radP[i]
		else:
		obj.power = 1.0

		def setCorePowerShape(me, radialPower=None, axialPower=None, coreStart=0.0):
		""" Use to set a power shape in the core region.

		@@ -1155,12 +1192,11 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		# renormalize
		p = list(np.array(p) / np.average(p))

		me.model.addAxialPowerShape('allrods', z, p)
		me.model.addAxialPowerShape('fuelrods', z, p)

		radPowers = []
		weights = []
		for obj in me.model.solidObjects.values():
		obj.powID = 'allrods'
		# Need to weight the rod power when normalizing radial power shape
		assert(obj.symtype == 1)
		if radialPower:
		@@ -1178,6 +1214,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		objType = me.model.solidTypeIDs[solidID]
		if isinstance(me.model.solidTypes[objType], Solid.FuelRod):
		obj.power = radialPower[i][j]
		obj.powID = 'fuelrods'
		else:
		raise RuntimeError(
		"Radial function not supported in setCorePowerShape")
		@@ -1186,15 +1223,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		radPowers.append(radialPower(
		utils.radialLocation(obj.x, obj.y)))

		# 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]
		me.normalizeRadialPowers()

		def setCorePowerShape3D(me, powerDist, wgts=None, coreStart=0.0, interp_type="nearest"):
		"""
		@@ -1242,15 +1271,7 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		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]
		me.normalizeRadialPowers()

		def _setCorePowerShapeAsm3D(me, powerDist, wgts=None, coreStart=0.0, interp_type="nearest"):
		"""

SubKit/build/genDeck.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -1124,7 +1124,7 @@ def writeDeck(model, filename):
		newDeckLines.append(makeGroupBanner("GROUP 17 - Channel/Rod Maps"))
		for l in group17Data:
		newDeckLines.append(l)
		if notrans == 1:
		if notrans == 1 or model.transPrintConv:
		newDeckLines.append(makeGroupBanner(
		"GROUP 19 - convergence Criteria for Steady State Solve"))
		for l in group19Data:

tests/nodal_p7/buildDeck.py

+54 −7

Original line number	Diff line number	Diff line
		@@ -10,6 +10,11 @@
		# - The ability to set percentTheoreticalDensity and the dynamic gap model is being tested
		# - The ability to set conductor options related to fuel material and gap models (Card 9.1 parameters)
		# - The ability to disable the main text output in Card Group 14
		# - The ability to provide the total core power to be used to calculate the average LHR
		# - The ability to set global multipliers on wide and narrow gap widths and gap length
		# - The ability to set RTWFP and convergence criteria for transient input decks
		# - The ability to set a core inlet temperature distribution
		# - The ability to set a core inlet flow rate distribution

		import sys
		sys.path.insert(0, '../../')
		@@ -30,13 +35,14 @@ def main():
		dGuideTubeInner = 0.561e-2 * 2 # m
		assemPitch = 21.5e-2 # m
		pitch = 1.26e-2 # m
		gapWidthMult = 0.5
		gapWidthMultWide = 0.1
		gapWidthMultNarrow = 0.5
		gapLengthMult = 2.0
		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
		corePow = 3.411e6 # kW
		directHeat = 0.026
		plenumPressure = 20.0 # bar
		plenumVolume = 8.783E-06 # m^3
		@@ -245,12 +251,16 @@ def main():
		builder = SquareLatticeLWR_Nodal(model, coreMap, pwr17x17, pitch, dz, solidGeoms, assemPitch=assemPitch,
		baffleGap=baffleGap, assemLosses=formLoss, assemLossLevels=gridLevels, rodDomains=rodDomains,
		modelGuideTubes=True,
		gapWidthMult=gapWidthMult, gapLengthMult=gapLengthMult)
		gapWidthMultWide=gapWidthMultWide, gapWidthMultNarrow=gapWidthMultNarrow, gapLengthMult=gapLengthMult)

		# Run to 1 second to let the solution reach steady state
		model.addTransientGroup(tEnd=1.0)
		# Test the ability to provide rtwfp during transients (used for system coupled calcs)
		model.addTransientGroup(tEnd=1.0, rtwfp=1e4)
		# Run a 10 second transient
		model.addTransientGroup(tEnd=11.0, editInterval=1.0)
		model.addTransientGroup(tEnd=11.0, editInterval=1.0, rtwfp=1e4)

		# Write card group 19 conv criteria to the input file, even for this transient model
		model.setSteadyState(transPrintConv=True, energyBalance=1e-4, massBalance=1e-5, maxIterations=200000)

		# Ramp the power down over time
		time = [1.0, 4.0, 4.3, 4.6, 4.9, 5.2, 5.5, 5.8, 6.1, 6.4, 6.7, 7.0, 11.0]
		@@ -264,11 +274,48 @@ def main():
		3.1, 3.4, 3.8, 4.0, 4.6, 5.4, 6.5, 7.9, 9.2, 10.0]
		flow = [1.00, 1.00, 0.90, 0.83, 0.78, 0.75, 0.71, 0.70, 0.67,
		0.66, 0.64, 0.62, 0.61, 0.58, 0.55, 0.51, 0.46, 0.43, 0.40]
		builder.setCoreInletBC(inletMassFlux, inletTemp, flowRamp=[time, flow])

		tempDist = [
		7, 6, 6, 7, 8, 8, 9,
		4, 3, 3, 3, 3, 4, 4, 6, 9, 10, 11,
		-2, 0, 0, -2, -1, 0, 0, 1, 5, 7, 9, 11, 12,
		-5, -5, -5, -5, -8, -6, -3, 0, 2, 7, 9, 11, 12,
		-9, -9, -11, -10, -12, -11, -10, -7, -1, 3, 6, 10, 11, 12, 13,
		-13, -14, -16, -16, -15, -11, -11, -8, -1, 3, 7, 10, 11, 12, 13,
		-17, -18, -18, -18, -14, -10, -6, -2, 0, 5, 6, 9, 11, 12, 13,
		-22, -22, -20, -17, -15, -8, -1, -1, -1, 5, 8, 10, 11, 12, 13,
		-26, -25, -21, -18, -12, -6, 0, 0, 4, 8, 9, 11, 12, 13, 13,
		-28, -26, -22, -15, -10, -3, 2, 5, 9, 9, 11, 11, 13, 13, 13,
		-27, -27, -20, -15, -6, -2, 3, 8, 9, 11, 11, 12, 13, 13, 13,
		-24, -20, -12, -5, -1, 4, 8, 10, 11, 12, 12, 12, 13,
		-18, -17, -9, -5, 0, 4, 7, 9, 10, 11, 11, 12, 12,
		-10, -7, -4, -1, 2, 4, 7, 9, 9, 10, 10,
		-3, -1, 1, 3, 5, 7, 7]

		flowDist = [
		0.85, 0.89, 0.88, 0.88, 0.86, 0.88, 0.80,
		0.80, 0.79, 0.89, 0.96, 1.04, 0.95, 0.90, 0.94, 0.81, 0.82, 0.76,
		0.79, 0.85, 0.80, 0.93, 1.11, 1.03, 1.08, 1.04, 1.14, 0.96, 0.73, 0.83, 0.82,
		0.83, 0.83, 0.98, 1.07, 1.07, 1.19, 1.06, 1.25, 1.10, 1.15, 1.00, 0.89, 0.87,
		0.82, 0.87, 0.95, 1.14, 1.08, 1.29, 1.11, 1.19, 1.09, 1.28, 1.09, 1.12, 0.99, 0.96, 0.85,
		0.83, 0.91, 1.03, 1.07, 1.32, 1.15, 1.28, 1.10, 1.19, 1.11, 1.20, 1.05, 1.05, 0.93, 0.87,
		0.85, 0.89, 0.99, 1.19, 1.12, 1.27, 1.13, 1.16, 1.09, 1.29, 1.07, 1.15, 1.00, 0.91, 0.85,
		0.85, 0.90, 0.88, 1.09, 1.27, 1.08, 1.15, 1.29, 1.12, 1.06, 1.15, 1.06, 0.86, 0.88, 0.84,
		0.85, 0.90, 0.95, 1.24, 1.10, 1.20, 1.07, 1.18, 1.12, 1.24, 1.10, 1.16, 0.93, 0.93, 0.87,
		0.83, 0.91, 0.98, 1.05, 1.25, 1.11, 1.22, 1.10, 1.26, 1.12, 1.26, 1.08, 1.06, 0.91, 0.85,
		0.83, 0.93, 0.95, 1.10, 1.09, 1.22, 1.12, 1.25, 1.11, 1.27, 1.08, 1.15, 0.95, 0.89, 0.85,
		0.86, 0.88, 0.99, 1.08, 1.08, 1.16, 1.09, 1.19, 1.06, 1.08, 0.94, 0.83, 0.83,
		0.81, 0.83, 0.86, 0.96, 1.03, 1.01, 1.05, 0.97, 0.98, 0.95, 0.78, 0.80, 0.80,
		0.80, 0.87, 0.90, 0.92, 0.98, 0.94, 0.97, 0.87, 0.85, 0.81, 0.79,
		0.85, 0.88, 0.91, 0.94, 0.90, 0.85, 0.82]


		builder.setCoreInletBC(inletMassFlux, inletTemp, flowRamp=[
		time, flow], massflowShape=flowDist, temperatureShape=tempDist)
		builder.setCoreOutletBC(outletPressure, inletTemp)

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

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

tests/nodal_p7/nodal_p7.inp.gold

+2497 −2456

File changed.

Preview size limit exceeded, changes collapsed.