Add ability to set axial variation tables (cff92c7b) · Commits · Gurecky, William / SubKit

SubKit/build/Model.py

+166 −2

Original line number	Diff line number	Diff line
		@@ -66,6 +66,20 @@ class Model:
		# which will either be the axial location or the power factor, depending on which dictionary it is.
		me.axialPowerShapeAxial = {}
		me.axialPowerShapePower = {}
		# Axial geometry modification tables. Of the form {ID: {'axial': [a1, a2, ...], 'mult': [m1, m2, ...]}}
		# Here ID is an ID for the table, the 'axial' list gives the axial level locations, and the 'mult' list
		# gives the multipliers for each level
		me.axialGeoShape = {}
		# Links the axial geometry modification tables to channels in the model. Of the form:
		# {chID: {'tables': [tab1, tab2, ...], 'aspects': [[aspect1], [aspect1, aspect2]]}}
		# A single channel may have multiple geometry modification tables applied to it, with each table applying
		# to multiple aspects of the geometry (momentum area, scalar area, wetted perimeter).
		me.channelGeoTables = {}
		# Links the axial geometry modification tables to gaps in the model. A gap can only have one
		# table applied, so the form is:
		# {gapID: tab1}
		# A geometry modification table can only apply to gap width.
		me.gapGeoTables = {}
		# List of FormLoss objects in the model
		me.formLosses = []
		# Same size as formLosses. Gives the channel ID of the loss
		@@ -229,6 +243,12 @@ class Model:
		# Under-relaxation coefficients
		me.betaHTC = None

		# Optional model specifications
		me.optionalModels = {}

		# Numerical controls
		me.numericalControls = {}

		@property
		def is_nodal(self):
		""" Access private is_nodal attr """
		@@ -252,7 +272,8 @@ class Model:
		"""
		me.title = str(title)

		def setModelOptions(me, beta=None, enableEntrainment=None, betaHTC=None):
		def setModelOptions(me, beta=None, enableEntrainment=None, betaHTC=None, intDragModel=None,
		flowRegimeMap=None, subcooledBoilingModel=None, tpFrictionModel=None):
		""" Set various modeling options

		Any model that can be changed via this method has a default.
		@@ -262,6 +283,15 @@ class Model:
		beta (float): Single-phase turbulent mixing coefficietn
		enableEntrainment (bool): Set to True to enable entrainment/de-entrainment of droplets
		betaHTC (float): Under-relaxation coefficient for fluid/solid heat transfer coefficient
		intDragModel (str): Enter "drift_flux" or "legacy" to pick the drift flux or legacy models for
		modeling interfacial drag, or leave as None to default.
		flowRegimeMap (str): Enter "legacy" or "ge_nonprop" to pick the legacy or GE non-proprietary
		model, or leave as None to default.
		subcooledBoilingModel (str): Enter "thom" to use Thom plus Hancox-Nicol for near-wall
		condensation (no ONB model) or "saha" to use the Gorenflo model with Saha-Zuber for determining
		bubble detachment (uses ONB model to determine when boiling starts).
		tpFrictionModel (str): Enter 'wallis', 'chisholm', or 'lockhart' to select from the Wallis, Chisholm
		or Lockhart-Martinelli models for the two-phase flow multiplier, or leave as None to default.
		"""

		if beta is not None:
		@@ -269,6 +299,18 @@ class Model:
		if enableEntrainment is not None:
		me.enableEntrainment = enableEntrainment
		me.betaHTC = betaHTC
		if intDragModel is not None:
		assert intDragModel in ['legacy', 'drift_flux']
		me.optionalModels['intDragModel'] = intDragModel
		if flowRegimeMap is not None:
		assert flowRegimeMap == 'legacy' or flowRegimeMap == 'ge_nonprop'
		me.optionalModels['flowRegimeMap'] = flowRegimeMap
		if subcooledBoilingModel is not None:
		assert subcooledBoilingModel == 'thom' or subcooledBoilingModel == 'saha'
		me.optionalModels['subcooledBoilingModel'] = subcooledBoilingModel
		if tpFrictionModel is not None:
		assert tpFrictionModel in ['wallis', 'chisholm', 'lockhart']
		me.optionalModels['tpFrictionModel'] = tpFrictionModel

		def setFrictionModel(me, model):
		""" Sets the friction model option
		@@ -1079,6 +1121,84 @@ class Model:
		me.axialPowerShapeAxial[time][powID] = axial
		me.axialPowerShapePower[time][powID] = power

		def addAxialGeoShape(me, ID, axialLevels, multipliers):
		"""
		Add a axially-dependent geometry modification table to the model

		This will specify a list of multipliers to be applied at axial levels specified in the
		axialLevels list. The axial levels shall be integer level indices. The meaning of the node
		(momentum level or scalar level) depends on how the table is applied. If the table is applied
		to momentum cell area, it will mean momentum level and if it is applied to scalar cell area, wetted
		perimeter, or gap width, it will mean scalar level.

		It is important to consider how CTF will apply these tables. A linear interpolation is done
		between points in the table to determine the local multiplier. For example, for the following table:

		1 0.5
		3 0.7
		4 0.8
		8 0.8

		The multiplier that will be applied to level 2 will be 0.6. However, if adjacent levels are
		specified, as in the case of 3 and 4 in the table, then no interpolation can be done. In this
		case, the 0.7 multiplier will be applied to Level 3 and 0.8 will be applied to Level 4. If
		the user desires no interpolation, then each region must be "capped", as in the case of Levels
		4 to 8. In this case, each level will have a value of 0.8 applied. Furthermore, the bottom and
		top of the model should be included in the table to ensure all multipliers are defined.

		Args:
		ID (int): Identifier for the table
		axialLevels (list): List of integer node levels
		multipliers (list): List of float multiplier values
		"""
		assert isinstance(ID, int)
		assert ID > 0
		assert len(axialLevels) == len(multipliers)
		assert all(multipliers) > 0.0
		assert ID not in me.axialGeoShape
		me.axialGeoShape[ID] = {'axial': axialLevels, 'mult': multipliers}

		def assignTableToChannel(me, tableID, chID, applies):
		""" Use this procedure to assign a geometry modification table to a channel in the model.

		The tableID must match a table ID that was defined using the :meth:`addAxialGeoTable` procedure.
		The chID must match a channel ID that was defined using the :meth:`addChannel` procedure.
		The `applies` input is a list of geometry aspects to which the table shall apply.
		The `applies` list can include "momentum", "scalar", and "pw", which represent the momentum
		cell areas, scalar cell areas, and wetted perimeter.

		Args:
		tableID (choice): ID of previously defined axial geometry table
		chID (int): ID of previously defined channel
		applies (list or str): Can provide a single string value or list of strings to define
		which aspects of the geometry this table shall apply
		"""
		if chID not in me.channelGeoTables:
		me.channelGeoTables[chID] = {'tables': [], 'aspects': []}
		me.channelGeoTables[chID]['tables'].append(tableID)
		validApplies = ['momentum', 'scalar', 'pw']
		if not isinstance(applies, list):
		ls = [applies]
		assert applies in validApplies
		else:
		ls = applies
		for val in applies:
		assert val in validApplies
		me.channelGeoTables[chID]['aspects'].append(ls)

		def assignTableToGap(me, tableID, gapID):
		""" Assigns a geometry modification table to a gap in the model

		The axial geometry modification table assigned to tableID and defined with :meth:`addAxialGeoTable`
		will be assigned to the specified gapID. The table multipliers will operate on the gap width.

		Args:
		tableID (choice): ID of previously defined geometry modification table
		gapID (int): ID of gap defined by :meth:`setGap` procedure
		"""
		assert gapID not in me.gapGeoTables
		me.gapGeoTables[gapID] = tableID

		def addTransientGroup(me, tEnd, editInterval=None, dtMin=None, dtMax=None):
		"""
		Adds a transient group to the simulation
		@@ -1132,6 +1252,18 @@ class Model:
		else:
		me.solver = 5

		def setNumericalControls(me, maxOuterIterations=None):
		""" Set numerical controls

		Args:
		maxOuterIterations (int): This is the maximum number of outer iterations to take during a timestep.
		Note that by default, CTF does not use the outer iteration loop. Adding this and setting number
		of outer iterations greater than 1 will enable the loop.
		"""
		if maxOuterIterations is not None:
		assert maxOuterIterations > 0
		me.numericalControls['maxOuterIterations'] = maxOuterIterations

		def setSteadyState(me, 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,
		@@ -1573,6 +1705,8 @@ class Model:

		me._validateBoundaryConditions()

		me._validateGeoTables()

		# Add a uniform power profile that covers the entire model so it can be assigned to rods
		# that had now power shape assigned
		if me._anyRodsWithoutPowerID():
		@@ -1772,7 +1906,6 @@ class Model:
		raise RuntimeError(
		"`setInitial` was not called and the channels do not all have inlet mass/temp and outlet pressure BCs with uniform values, so cannot infer the correct initial conditions")


		# Determine the start/end axial locations of the sections
		z = 0.0
		for secID in sorted(me.sections.keys()):
		@@ -2088,6 +2221,37 @@ class Model:
		raise RuntimeError("For solidID: " + str(solidID) + " in Section: " + str(
		u) + " percentages of connected channels add up to " + str(totPercent) + " instead of 1.0")

		def _validateGeoTables(me):
		""" Ensure all geometry tables applied to model objects were defined and that they were all assigned
		to objects that exist in the model.
		"""
		# Ensure axial geo table IDs are sequential
		ID = 1
		for key in sorted(me.axialGeoShape.keys()):
		if key != ID:
		raise RuntimeError(
		"Axial geometry table IDs must be sequential starting at 1.")
		ID = ID + 1
		for chID in me.channelGeoTables.keys():
		foundChannel = False
		for section in me.sections.keys():
		if chID in me.sections[section].channels.keys():
		foundChannel = True
		if not foundChannel:
		raise RuntimeError("Geometry table ID: " + str(me.channelGeoTables[chID]['tables'][0]) +
		" assigned to non-existent channel: " + str(chID))
		for tableID in me.channelGeoTables[chID]['tables']:
		if tableID not in me.axialGeoShape:
		raise RuntimeError("Geometry table ID: " + str(tableID) + " was not defined, but is trying "
		+ "to be applied to channel: " + str(chID))
		for gapID in me.gapGeoTables.keys():
		if gapID not in me.gaps:
		raise RuntimeError("Geometry table ID: " + str(me.gapGeoTables[gapID]) + " attempting to be applied"
		+ " to non-existent gap: " + str(gapID))
		if tableID not in me.axialGeoShape:
		raise RuntimeError("Geometry table ID: " + str(tableID) + " was not defined, but is trying "
		+ "to be applied to gap: " + str(gapID))


		class Gap(object):
		""" Defines a gap (a lateral connection between two channels in the same axial section)

SubKit/build/genDeck.py

+92 −3

Original line number	Diff line number	Diff line
		@@ -85,7 +85,10 @@ def writeDeck(model, filename):
		controlData.append("*INITIAL DUMPF\n")
		controlData.append(" 1 0\n")
		controlData.append("** EPSO OITMAX IITMAX COURANT\n")
		controlData.append(" 1e-4 5 200 0.800000\n")
		oitmax = 5
		if 'maxOuterIterations' in model.numericalControls:
		oitmax = model.numericalControls['maxOuterIterations']
		controlData.append(" 1e-4{:14d} 200 0.800000\n".format(oitmax))
		controlData.append("*TITLE\n")
		controlData.append(model.title + "\n")

		@@ -158,6 +161,21 @@ def writeDeck(model, filename):
		if model.betaHTC is not None:
		group1Data.append("{{beta_htc}} {:15.3e}\n".format(model.betaHTC))

		if 'intDragModel' in model.optionalModels:
		group1Data.append("{{int_drag_model}} {:s}\n".format(
		model.optionalModels['intDragModel']))
		if 'flowRegimeMap' in model.optionalModels:
		group1Data.append("{{flow_regime_map}} {:s}\n".format(
		model.optionalModels['flowRegimeMap']))
		if 'subcooledBoilingModel' in model.optionalModels:
		if model.optionalModels['subcooledBoilingModel'] == 'saha':
		group1Data.append("{use_onb_model} .true.\n")
		if 'tpFrictionModel' in model.optionalModels:
		group1Data.append("{{tp_fric_model}} {:s}\n".format(
		model.optionalModels['tpFrictionModel']))
		if 'maxOuterIterations' in model.numericalControls:
		group1Data.append("{enable_outer_iteration_loop} .true.\n")

		# The Card Group 2 data to be written to the deck
		group2Data = []
		group2Data.append("**NGR\n")
		@@ -295,6 +313,63 @@ def writeDeck(model, filename):
		group4Data.append("** MSIM\n")
		group4Data.append("{0:8d}\n".format(model._getNCELLS()))

		group5Data = []
		group6Data = []
		if len(model.gapGeoTables) > 0 or len(model.channelGeoTables) > 0:
		group5Data.append("**NGR\n5\n")
		group5Data.append("*Card 5.1\n")
		maxLen = 0
		for key in model.axialGeoShape.keys():
		maxLen = max(maxLen, len(model.axialGeoShape[key]['axial']))
		group5Data.append("** NAFACT NAXL\n")
		group5Data.append("{:9d}{:7d} 0 0 0 0 0 0 0 0 0 0 0 0\n".format(
		len(model.axialGeoShape.keys()), maxLen))
		group5Data.append("**Card 5.2\n")
		for key in sorted(model.axialGeoShape.keys()):
		line = ""
		for i in range(len(model.axialGeoShape[key]['axial'])):
		line = line + "{:5d}{:16.4e}".format(model.axialGeoShape[key]['axial'][i],
		model.axialGeoShape[key]['mult'][i])
		# Pad the rest with zeros if this table is shorter than max
		for i in range(len(model.axialGeoShape[key]['axial']), maxLen):
		line = line + "{:5d}{:16.4e}".format(0, 0.0)
		group5Data.append(line + "\n")

		group6Data.append("**NGR\n6\n")
		group6Data.append("*Card 6.1\n")
		group6Data.append("** N1\n")
		group6Data.append("{:7d} 0 0 0 0 0 0 0 0 0 0 0 0 0\n".format(len(model.channelGeoTables.keys())
		+ len(model.gapGeoTables.keys())))
		group6Data.append("*Card 6.2\n")
		group6Data.append("** IACT IAMT IPWT ICRG\n")

		def getMom(channelGeoTables):
		""" Takes the channelGeoTables entry for chID and returns the flag for IAMT in the CTF input deck."""
		for i, tableID in enumerate(channelGeoTables['tables']):
		if 'momentum' in channelGeoTables['aspects'][i]:
		return tableID
		return 0

		def getScal(channelGeoTables):
		""" Takes the channelGeoTables entry for chID and returns the flag for IACT in the CTF input deck."""
		for i, tableID in enumerate(channelGeoTables['tables']):
		if 'scalar' in channelGeoTables['aspects'][i]:
		return tableID
		return 0

		def getPw(channelGeoTables):
		""" Takes the channelGeoTables entry for chID and returns the flag for IPWT in the CTF input deck."""
		for i, tableID in enumerate(channelGeoTables['tables']):
		if 'pw' in channelGeoTables['aspects'][i]:
		return tableID
		return 0
		for chID in sorted(model.channelGeoTables.keys()):
		group6Data.append("{:7d}{:7d}{:7d}{:7d} 0 0 0 0 0 0 0 0 0 0 0\n".format(getScal(model.channelGeoTables[chID]),
		getMom(model.channelGeoTables[chID]), getPw(model.channelGeoTables[chID]), chID))
		for gapID in sorted(model.gapGeoTables.keys()):
		group6Data.append("{:7d}{:7d}{:7d}{:7d} 0 0 0 0 0 0 0 0 0 0 0\n".format(-model.gapGeoTables[gapID], 0, 0,
		gapID))

		group7Data = []
		group7Data.append("**NGR\n")
		group7Data.append(" 7\n")
		@@ -409,8 +484,12 @@ def writeDeck(model, filename):
		dnbchk = 1
		elif model.conductorOptions['postChfCheck']['value'] == 'groeneveld':
		dnbchk = 3
		group8Data.append("{:7d}{:7d}{:6d} 1 0 0 1 1 0 0{:6d} 1{:6d} 0\n".format(
		len(heated), len(unheated), model.conductorOptions['nc']['value'], w3, dnbchk))
		ihtc = 1
		if 'subcooledBoilingModel' in model.optionalModels:
		if model.optionalModels['subcooledBoilingModel'] == 'saha':
		ihtc = 2
		group8Data.append("{:7d}{:7d}{:6d} 1 0 0 1 1 0 0{:6d}{:5d}{:6d} 0\n".format(
		len(heated), len(unheated), model.conductorOptions['nc']['value'], w3, ihtc, dnbchk))

		# If this is a parallel model write data
		if model.rodsInDomain:
		@@ -995,6 +1074,16 @@ def writeDeck(model, filename):
		"GROUP 4 - Vertical Channel Connection Data"))
		for l in group4Data:
		newDeckLines.append(l)
		if len(group5Data) > 0:
		newDeckLines.append(makeGroupBanner(
		"GROUP 5 - Axial Geometry Variation Tables"))
		for l in group5Data:
		newDeckLines.append(l)
		if len(group6Data) > 0:
		newDeckLines.append(makeGroupBanner(
		"GROUP 6 - Channel/Gap Variation Table Application"))
		for l in group6Data:
		newDeckLines.append(l)
		newDeckLines.append(makeGroupBanner(
		"GROUP 7 - Grid Loss Coefficient Data"))
		for l in group7Data:

doc/input_file.rst

+2 −2

Original line number	Diff line number	Diff line
		@@ -448,13 +448,13 @@ temp

		enth
		----
		:Description: Initial enthalpy for the system
		:Description: Initial enthalpy for the system - Feature to be implemented in a future release.
		:Min Occurance: 0
		:Max Occurance: 1
		:Units: kJ/kg
		:Type: Real
		:Range: >0.0
		:Note: If enthalpy is provided, temperature must still be specified for setting the initial temperature in solid objects
		:Note: If enthalpy is provided along with temperature, the enthalpy will apply to the liquid and temperature will be used for the solids.

		qprime
		------

tests/inputRegression/vuq_mult.txt

0 → 100644

+8 −0

Original line number	Diff line number	Diff line
		k_hspl = 1.000000e+00
		k_xk_anrflm = 1.000000e+00
		k_xk_sb = 1.000000e+00
		k_xk_slg = 1.000000e+00
		ka_hspl = 0.000000e+00
		ka_xk_anrflm = 0.000000e+00
		ka_xk_sb = 0.000000e+00
		ka_xk_slg = 0.000000e+00

tests/inputRegression/vuq_param.txt

0 → 100644

+13 −0

Original line number	Diff line number	Diff line
		k_db_1 = 2.300000e-02
		k_db_2 = 8.000000e-01
		k_db_3 = 4.000000e-01
		k_db_4 = 7.860000e+00
		k_tke_coeff = 3.000000e-01
		k_xkwlxA = 1.100000e+00
		k_xkwlxB = -2.200000e+00
		k_xkwlxChA = 1.000000e+00
		k_xkwlxChB = 1.000000e+00
		k_xkwlxChC = 2.000000e+00
		k_xkwlxChD = 1.000000e+00
		k_xkwlxChE = 1.000000e+00
		k_xkwlxChF = 1.000000e+10