Add support for specifying enthalpy as IC and BC (73a05f17) · Commits · Gurecky, William / SubKit

SubKit/build/InpBuild.py

+9 −5

Original line number	Diff line number	Diff line
		@@ -47,12 +47,16 @@ class InpBuilder(object):
		mdotinit = me.inp.getInitialmdot()
		pressureinit = me.inp.getInitialPressure()
		temperatureinit = me.inp.getInitialTemperature()
		#enthalpyinit = me.inp.getInitialEnthalpy()
		enthalpyinit = me.inp.getInitialEnthalpy()

		if massfluxinit:
		mdotinit = massfluxinit * me.inp.getInFlowArea()
		me.model.setInitialConditions(
		mdotInit=mdotinit, tempInit=temperatureinit, pressureInit=pressureinit)
		if enthalpyinit is not None:
		temp_solid = temperatureinit
		temperatureinit = None
		else:
		temp_solid = None

		me.model.setInitial(mdot=mdotinit, massflux=massfluxinit, temp=temperatureinit,
		temp_solid=temp_solid, enthalpy=enthalpyinit, pressure=pressureinit)

		# Boundary Conditions
		for chID in me.inp.getChannelIDs():

SubKit/build/InpParse.py

+13 −3

Original line number	Diff line number	Diff line
		@@ -195,6 +195,10 @@ class InpParse(object):
		if 'temp' in me.inpDict['initial']:
		return me.inpDict['initial']['temp']['value']

		def getInitialEnthalpy(me):
		if 'enth' in me.inpDict['initial']:
		return me.inpDict['initial']['enth']['value']

		def getSectionIDs(me):
		return me.secIDs

		@@ -414,11 +418,17 @@ class InpParse(object):
		massflux = me.getInitialMassFlux()
		mdotin = me.getInitialmdot()
		temperature = me.getInitialTemperature()
		enthalpy = me.getInitialEnthalpy()
		pressure = me.getInitialPressure()
		if massflux:
		mdot = massflux * me.getChannelArea(chID)
		if mdotin:
		mdot = mdotin * me.getChannelArea(chID) / me.getInFlowArea()
		if enthalpy:
		BCs.append(BC.MassEnthalpy(mdot=mdot, h=enthalpy))
		BCs.append(BC.PressureEnthalpy(
		pressure=pressure, h=enthalpy))
		else:
		BCs.append(BC.MassTemperature(mdot=mdot, T=temperature))
		BCs.append(BC.PressureTemperature(
		pressure=pressure, T=temperature))

SubKit/build/Model.py

+98 −68

Original line number	Diff line number	Diff line
		@@ -424,39 +424,67 @@ class Model:
		mdot (float): Initial mass flow rate in the system
		massflux (float): Initial mass flux in the system
		temp (float): Initial temperature in the model
		temp_solid (float): Initial temperature of solids in the model
		enthalpy (float): Initial enthalpy in the model
		pressure (float): Initial (and reference) pressure in the system

		"""
		if me.setInitialCalled:
		raise RuntimeError("`setInitial` already called")
		assert not ('mdot' in kwargs and 'massflux' in kwargs)
		assert 'mdot' in kwargs or 'massflux' in kwargs
		assert not ('temp' in kwargs and 'enthalpy' in kwargs)
		assert 'temp' in kwargs or 'enthalpy' in kwargs
		if 'enthalpy' in kwargs:
		assert 'temp_solid' in kwargs
		assert 'pressure' in kwargs

		if 'mdot' in kwargs:
		me.mdotInit = kwargs['mdot']
		mdot = kwargs['mdot']
		else:
		mdot = None
		if 'massflux' in kwargs:
		massflux = kwargs['massflux']
		else:
		massflux = None
		if 'temp' in kwargs:
		temp = kwargs['temp']
		else:
		temp = None
		if 'temp_solid' in kwargs:
		temp_solid = kwargs['temp_solid']
		else:
		temp_solid = None
		if 'enthalpy' in kwargs:
		enthalpy = kwargs['enthalpy']
		else:
		enthalpy = None
		if 'pressure' in kwargs:
		pressure = kwargs['pressure']
		else:
		pressure = None

		# Cannot have both mdot and massflux provided
		assert not ((mdot is not None) and (massflux is not None))
		# Cannot have both temp and enthalpy provided
		assert not ((temp is not None) and (enthalpy is not None))
		# Must specify at least temperature or enthalpy
		assert (temp is not None) or (enthalpy is not None)
		if enthalpy is not None:
		assert temp_solid is not None
		assert pressure is not None

		if mdot is not None:
		me.mdotInit = mdot
		me.massFluxInit = me.mdotInit / me.sections[1].getFlowArea()
		me.mdotInitialWasSpecified = True
		else:
		me.massFluxInit = kwargs['massflux']
		me.massFluxInit = massflux
		me.mdotInit = me.massFluxInit * me.sections[1].getFlowArea()

		if 'temp' in kwargs:
		me.tempInit = kwargs['temp']
		if temp is not None:
		me.tempInit = temp
		else:
		me.enthalpyInit = kwargs['enthalpy']
		me.enthalpyInit = enthalpy

		if 'temp_solid' in kwargs:
		me.solidTempInit = kwargs['temp_solid']
		if temp_solid is not None:
		me.solidTempInit = temp_solid
		else:
		me.solidTempInit = kwargs['temp']
		me.solidTempInit = temp

		me.pressureInit = kwargs['pressure']
		me.pressureInit = pressure

		me.setInitialCalled = True

		@@ -925,12 +953,14 @@ class Model:
		if all([isinstance(entry, list) for entry in gridEnhancementCoeffs]):
		assert(len(gridEnhancementCoeffs) == len(levelsList))
		for coeffs in gridEnhancementCoeffs:
		assert(all([isinstance(coeff,float) for coeff in coeffs]))
		assert(all([isinstance(coeff, float)
		for coeff in coeffs]))
		assert(len(coeffs) == 4)
		gridEnhancementCoeffsList = gridEnhancementCoeffs
		elif all([isinstance(entry, float) for entry in gridEnhancementCoeffs]):
		assert(len(gridEnhancementCoeffs) == 4)
		gridEnhancementCoeffsList = [gridEnhancementCoeffs] * len(levelsList)
		gridEnhancementCoeffsList = [
		gridEnhancementCoeffs] * len(levelsList)
		else:
		raise RuntimeError(
		"gridEnhancementCoeffs must be a list of lists or a list of floats " + str(gridEnhancementCoeffs))
		@@ -1316,12 +1346,12 @@ class Model:
		if type(value) is not type(me.conductorOptions[name]['value']):
		# Assume all entries must have the same type as the
		# default entry set in this class.
		raise TypeError(name+" must be "+str(type(me.conductorOptions[name]['value']))+
		" but it was provided as "+str(type(value)))
		raise TypeError(name + " must be " + str(type(me.conductorOptions[name]['value']))
		+ " but it was provided as " + str(type(value)))
		if value not in me.conductorOptions[name]['validEntries']:
		raise ValueError("The value provided for " + name + " (" + str(value) +
		") is not one of the valid values: "
		+str(me.conductorOptions[name]['validEntries']))
		") is not one of the valid values: " +
		str(me.conductorOptions[name]['validEntries']))
		me.conductorOptions[name]['value'] = value
		me.conductorOptions[name]['set'] = True

SubKit/build/SquareLatticeLWR_Nodal.py

+29 −24

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


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

		@@ -796,8 +797,8 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):

		# Convert the passed assemLossLevels to a dictionary that takes assembly type ID and returns
		# the levels in that location.
		assert(isinstance(assemLossLevels, list)
		or isinstance(assemLossLevels, dict))
		assert(isinstance(assemLossLevels, list) or
		isinstance(assemLossLevels, dict))
		if isinstance(assemLossLevels, list):
		assemLossLevelsDict = {}
		for assemTypeIdx in chMapObj.getAssemTypeIds():
		@@ -808,16 +809,19 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		assert(assemTypeIdx in assemLossLevels)
		assemLossLevelsDict = assemLossLevels

		assemLossesDict = me.convertAssemInputToDict(assemLosses,assemLossLevelsDict,chMapObj)
		assemLossesDict = me.convertAssemInputToDict(
		assemLosses, assemLossLevelsDict, chMapObj)
		# Grid blockage and enhancement are optional if loss levels were provided.
		if assemGridBlockageRatios is not None:
		assemGridBlockageRatiosDict = me.convertAssemInputToDict(assemGridBlockageRatios,assemLossLevelsDict,chMapObj)
		assemGridBlockageRatiosDict = me.convertAssemInputToDict(
		assemGridBlockageRatios, assemLossLevelsDict, chMapObj)
		else:
		assemGridBlockageRatiosDict = {}
		for assemTypeIdx in chMapObj.getAssemTypeIds():
		assemGridBlockageRatiosDict[assemTypeIdx] = None
		if assemGridEnhancementCoeffs is not None:
		assemGridEnhancementCoeffsDict = me.convertAssemListInputToDict(assemGridEnhancementCoeffs,assemLossLevelsDict,chMapObj)
		assemGridEnhancementCoeffsDict = me.convertAssemListInputToDict(
		assemGridEnhancementCoeffs, assemLossLevelsDict, chMapObj)
		else:
		assemGridEnhancementCoeffsDict = {}
		for assemTypeIdx in chMapObj.getAssemTypeIds():
		@@ -846,7 +850,8 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):

		# Grid blockages and enhancement coeffs must be defined
		# for the same assembly type IDs.
		assert(assemGridBlockageRatiosDict.keys()==assemGridEnhancementCoeffsDict.keys())
		assert(assemGridBlockageRatiosDict.keys()
		== assemGridEnhancementCoeffsDict.keys())
		for assemTypeID in assemGridBlockageRatiosDict.keys():
		# Find all chans associated with this assembly type
		# and assign the grid enhancement parameters for them.
		@@ -1004,8 +1009,8 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		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")
		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)

		@@ -1028,8 +1033,8 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		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.")
		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)
		@@ -1316,8 +1321,8 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		asm_axial_power_profile[z_lvl]

		# Normalize such that axially integrated asm nodal power is 1
		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))
		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))

		# set axial power profile in the ctf model
		me.model.addAxialPowerShape(
		@@ -1433,8 +1438,8 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
		node_axial_power_profile[z_lvl]

		# Normalize such that axially integrated asm nodal power is 1
		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))
		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(

SubKit/build/genDeck.py

+21 −16

Original line number	Diff line number	Diff line
		@@ -170,8 +170,8 @@ def writeDeck(model, filename):

		# If this is a parallel model
		if model.channelsInDomain:
		group2Data.append("{num_domain} " +
		str(len(model.channelsInDomain)) + "\n")
		group2Data.append("{num_domain} "
		+ str(len(model.channelsInDomain)) + "\n")
		for domain in sorted(model.channelsInDomain.keys()):
		group2Data.append("{domain} " + str(domain) + "\n")
		group2Data.append("** local_inex global_index\n")
		@@ -315,7 +315,8 @@ def writeDeck(model, filename):
		len(model.formLosses), model.numGridTypes, ifespv))
		if model.numGridTypes > 0:
		group7Data.append("*Card 7.1.5\n")
		group7Data.append("** C1 C2 C3 C4 A PHI ABLOC SPBLOC\n")
		group7Data.append(
		"** C1 C2 C3 C4 A PHI ABLOC SPBLOC\n")
		for grid in range(model.numGridTypes):
		# Use default values; overwrite if values for this grid are stored in model.gridEnhancementDict
		gridCoeffs = [5.55, -0.13, -0.034, 0.4]
		@@ -332,7 +333,8 @@ def writeDeck(model, filename):
		gridCoeffs[0], gridCoeffs[1], gridCoeffs[2], gridCoeffs[3], abloc))
		if model.numGridTypes > 0:
		group7Data.append("*Card 7.1.6\n")
		group7Data.append("** GRD J CD1 CD2 CD3 CD4 CD5 CD6 CD7 CD8 CD9 CD10 CD11 CD12 \n")
		group7Data.append(
		"** GRD J CD1 CD2 CD3 CD4 CD5 CD6 CD7 CD8 CD9 CD10 CD11 CD12 \n")
		for grid in range(model.numGridTypes):
		chanList = model.gridEnhancementDict[grid]['chans']
		ichan = -1
		@@ -412,8 +414,8 @@ def writeDeck(model, filename):

		# If this is a parallel model write data
		if model.rodsInDomain:
		group8Data.append("{num_domain} " +
		str(len(model.rodsInDomain)) + "\n")
		group8Data.append("{num_domain} "
		+ str(len(model.rodsInDomain)) + "\n")
		for domain in sorted(model.rodsInDomain.keys()):
		group8Data.append("{domain} " + str(domain) + "\n")
		group8Data.append("** local_inex global_index owner\n")
		@@ -559,8 +561,10 @@ def writeDeck(model, filename):
		numtype = len(list(set(model.solidTypeIDs.values())))
		group9Data.append(
		"{0:7d}{1:5d}{2:5d}{3:5d}{4:5d}{5:5d}{6:5d}{7:5d} 0 0 0 0 0 0\n".format(
		numtype, model.conductorOptions['irelf']['value'], model.conductorOptions['iconf']['value'], model.conductorOptions['imwr']['value'],
		model.conductorOptions['ifswell']['value'], model.conductorOptions['ifdens']['value'], model.conductorOptions['ifreloc']['value'],
		numtype, model.conductorOptions['irelf']['value'], model.conductorOptions[
		'iconf']['value'], model.conductorOptions['imwr']['value'],
		model.conductorOptions['ifswell']['value'], model.conductorOptions[
		'ifdens']['value'], model.conductorOptions['ifreloc']['value'],
		model.conductorOptions['ifcladcreep']['value']))

		def charID(pinobj):
		@@ -591,7 +595,8 @@ def writeDeck(model, filename):
		if igpc == -1:
		# Default VERAIn values for Problem 7
		group9Data.append('*Card 9.3\n')
		group9Data.append('** PGAS VPLEN ROUFF ROUFC HE XE AR KR H N OXIDET\n')
		group9Data.append(
		'** PGAS VPLEN ROUFF ROUFC HE XE AR KR H N OXIDET\n')
		group9Data.append('{:12.3e}{:12.3e} -1.000E-02 -1.000E-02 1.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00\n'.format(
		typeobj.plenumPressure, typeobj.plenumVolume))
		group9Data.append('*Card 9.4\n')