Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow components to have children in support of particle fuel modeling #702

Open
wants to merge 7 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
61 changes: 49 additions & 12 deletions armi/reactor/blocks.py
Original file line number Diff line number Diff line change
Expand Up @@ -276,7 +276,7 @@ def getSmearDensity(self, cold=True):
# Compute component areas
cladID = numpy.mean([clad.getDimension("id", cold=cold) for clad in clads])
innerCladdingArea = (
math.pi * (cladID ** 2) / 4.0 * self.getNumComponents(Flags.FUEL)
math.pi * (cladID**2) / 4.0 * self.getNumComponents(Flags.FUEL)
)
fuelComponentArea = 0.0
unmovableComponentArea = 0.0
Expand Down Expand Up @@ -751,6 +751,43 @@ def _updateDetailedNdens(self, frac, adjustList):
# on the interface stack.
self.p.pdensDecay *= frac

def updateComponentGroupMults(self):
"""
Update component group mults.

This applies only to block child components that themselves have child components.

They were temporarily set to the desired blend fraction during component
construction. This derives the actual multiplicity based on the target
blend fraction.

Note that the blend fractions are applied to the hot thermally-expanded
dimensions.

In order to account properly for the mass of the background material
in addition to the materials of the children, we must also shrink
the background component's mult.

See Also
--------
reactor.blueprints.componentBlueprint.ComponentBlueprint : sets the mults
reactor.blurprints.blockBlueprint.BlockBlueprint.construct : calls this
"""
for c in self.getChildren():
if len(c) > 1:
backgroundVol = c.getVolume()
# expect all mults to be temporarily set to blend frac from blueprints
blendFrac = c[0].p.mult
# remove scaling by mult=blendFrac here
childVol = sum(subchild.getVolume() / blendFrac for subchild in c)
# mult should be such that childVol/backgroundVol = blendFrac
# so: childVol * newMult = blendFrac * backgroundVol
newMult = blendFrac * backgroundVol / childVol
for subchild in c:
subchild.setDimension("mult", newMult)
# shrink background component to fit
c.setDimension("mult", c.p.mult * (1 - blendFrac))

def completeInitialLoading(self, bolBlock=None):
"""
Does some BOL bookkeeping to track things like BOL HM density for burnup tracking.
Expand Down Expand Up @@ -1905,14 +1942,10 @@ def getSymmetryFactor(self):
# seeing the first one is the easiest way to detect them.
# Check it last in the and statement so we don't waste time doing it.
upperEdgeLoc = self.core.spatialGrid[-1, 2, 0]
if (
symmetryLine
in [
grids.BOUNDARY_0_DEGREES,
grids.BOUNDARY_120_DEGREES,
]
and bool(self.core.childrenByLocator.get(upperEdgeLoc))
):
if symmetryLine in [
grids.BOUNDARY_0_DEGREES,
grids.BOUNDARY_120_DEGREES,
] and bool(self.core.childrenByLocator.get(upperEdgeLoc)):
return 2.0
return 1.0

Expand Down Expand Up @@ -1950,11 +1983,15 @@ def autoCreateSpatialGrids(self):

ringNumber = hexagon.numRingsToHoldNumCells(self.getNumPins())
# For the below to work, there must not be multiple wire or multiple clad types.
pitch = self.getPinPitch(cold=True)
if pitch is None:
raise ValueError(
f"Could not create spatialGrid for block {self.p.type} "
"because the pin pitch was not detected."
)
# note that it's the pointed end of the cell hexes that are up (but the
# macro shape of the pins forms a hex with a flat top fitting in the assembly)
grid = grids.HexGrid.fromPitch(
self.getPinPitch(cold=True), numRings=0, pointedEndUp=True
)
grid = grids.HexGrid.fromPitch(pitch, numRings=0, pointedEndUp=True)
spatialLocators = grids.MultiIndexLocation(grid=self.spatialGrid)
numLocations = 0
for ring in range(ringNumber):
Expand Down
6 changes: 6 additions & 0 deletions armi/reactor/blueprints/blockBlueprint.py
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,7 @@ def construct(
Top layer groups the by-block material modifications under the `byBlock` key
and the by-component material modifications under the component's name.
The inner dict under each key contains material modification names and values.

"""
runLog.debug("Constructing block {}".format(self.name))
components = collections.OrderedDict()
Expand Down Expand Up @@ -194,6 +195,11 @@ def construct(
b.autoCreateSpatialGrids()
except (ValueError, NotImplementedError) as e:
runLog.warning(str(e), single=True)

# now that components are in blocks we have heights and can actually
# compute the mults of the component groups.
b.updateComponentGroupMults()

return b

def _checkByComponentMaterialInput(self, materialInput):
Expand Down
48 changes: 45 additions & 3 deletions armi/reactor/blueprints/componentBlueprint.py
Original file line number Diff line number Diff line change
Expand Up @@ -153,6 +153,8 @@ def shape(self, shape): # pylint: disable=no-self-use; reason=yamlize requireme
isotopics = yamlize.Attribute(type=str, default=None)
latticeIDs = yamlize.Attribute(type=list, default=None)
origin = yamlize.Attribute(type=list, default=None)
blends = yamlize.Attribute(type=list, default=None)
blendFracs = yamlize.Attribute(type=list, default=None)
orientation = yamlize.Attribute(type=str, default=None)
mergeWith = yamlize.Attribute(type=str, default=None)
area = yamlize.Attribute(type=float, default=None)
Expand All @@ -175,9 +177,49 @@ def construct(self, blueprint, matMods):
constructedObject.add(component)

else:
constructedObject = components.factory(shape, [], kwargs)
_setComponentFlags(constructedObject, self.flags, blueprint)
insertDepletableNuclideKeys(constructedObject, blueprint)
if self.blends:
# build a blended object
for groupName, blendFrac in zip(self.blends, self.blendFracs):
# blendFrac is a volume fraction, and so we need to adjust the multiplicities
# so that the background component and the child group match up
# The area of the background component will be fully determined by its dims and mult.
# but internally, its number densities and volumes need to be set to match
# the blendFrac
group = blueprint.componentGroups[groupName]
# strip off the blend/blendFrac args since they aren't valid on any
# specific shape's constructor
del kwargs["blends"]
del kwargs["blendFracs"]
# build the background object
constructedObject = components.factory(shape, [], kwargs)
# build the child objects
# all grouped components have to be input with the same mult for now
inputMult = None
for groupedComponent in group:
componentDesign = blueprint.componentDesigns[
groupedComponent.name
]
if inputMult is None:
inputMult = componentDesign.mult
if componentDesign.mult != inputMult:
raise ValueError(
f"Grouped components must all have the same input mult of {inputMult}"
)
component = componentDesign.construct(blueprint, matMods=dict())
# temporarily set grouped component mults to the blend fraction
# these will be updated based on parent component volume
# during block construction
component.setDimension("mult", blendFrac)
_setComponentFlags(component, self.flags, blueprint)
insertDepletableNuclideKeys(component, blueprint)
constructedObject.add(component)
children = {c.name: c for c in constructedObject.getChildren()}
for child in children.values():
child.resolveLinkedDims(children)
Comment on lines +216 to +218
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This enables individual layers to have linked dimensions, right? Nice!

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yeah it should. though we should make more unit tests proving it.

else:
constructedObject = components.factory(shape, [], kwargs)
_setComponentFlags(constructedObject, self.flags, blueprint)
insertDepletableNuclideKeys(constructedObject, blueprint)
return constructedObject

def _conformKwargs(self, blueprint, matMods):
Expand Down
64 changes: 44 additions & 20 deletions armi/reactor/components/component.py
Original file line number Diff line number Diff line change
Expand Up @@ -609,30 +609,42 @@ def getNuclides(self):

This includes anything that has been specified in here, including trace nuclides.
"""
return list(self.p.numberDensities.keys())
nucs = set(self.p.numberDensities.keys())
if self._children:
childNuclides = set(composites.Composite.getNuclides(self))
nucs |= childNuclides
return nucs
Comment on lines -612 to +616
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same concern on getting nuclides just in the component vs. component and all the children

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yeah that's a good concern. I'd hate to make you basically re-write these methods for accessing the details of the matrix itself. We could add in a flag or something as a method argument? getNuclides(includeChildren=False) or something?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

includeChildren=False feels like a good solution. Then whatever application or plugin is doing the writing could check if a Component has children prior to writing the material definitions


def getNumberDensity(self, nucName):
def getNuclideNumberDensities(self, nucNames):
"""
Get the number density of nucName, return zero if it does not exist here.
Return a list of number densities for the nuc names requested.

Parameters
----------
nucName : str
Nuclide name
If this Component has children, then homogenize their number densities in.

Returns
-------
number density : float
number density in atoms/bn-cm.
Component ndens is unique in that it combines its own actual composition
with that of its potential children
"""
return self.p.numberDensities.get(nucName, 0.0)

def getNuclideNumberDensities(self, nucNames):
"""Return a list of number densities for the nuc names requested."""
return [self.p.numberDensities.get(nucName, 0.0) for nucName in nucNames]

def _getNdensHelper(self):
return dict(self.p.numberDensities)
if self._children:
childDens = dict(
zip(
nucNames,
composites.Composite.getNuclideNumberDensities(self, nucNames),
)
)
# get volume of children by using the composite method which loops over all children
childVol = composites.Composite.getVolume(self)
# get self volume by calling the shape-specific volume method of the background shape
totalVol = self.getVolume()
childFrac = childVol / totalVol
selfFrac = 1.0 - childFrac
return [
self.p.numberDensities.get(nucName, 0.0) * selfFrac
+ childDens.get(nucName, 0.0) * childFrac
for nucName in nucNames
]
else:
return [self.p.numberDensities.get(nucName, 0.0) for nucName in nucNames]
Comment on lines +628 to +647
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How would a plugin be able to get the number densities of just the matrix component with this change?

On one hand this is consistent with how the design of the composite tree: Composite.getNuclideNumberDensities gives the densities for all children on a node. So this would match well with Block.getNuclideNumberDensities for some nodal homogenization methods. But modeling this Component exactly in a monte carlo code would be problematic w/o the number densities of just the Component and not its children

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great point. You would have to access the matrix's .p.numberDensities parameter directly in order to differentiate between it and its children with this change as written. Willing to consider more elegant options.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ntouran is

matrix's .p.numberDensities

still the recommended approach here? I've been testing this out over the last few days and had a hard time remembering this thread.

Would something like Component.getNuclideNumberDensities(children=False) be acceptable? I tend to be hesitant pulling data from these parameter containers as the feel (maybe incorrectly) more private than a public method


def setName(self, name):
"""Components use name for type and name."""
Expand Down Expand Up @@ -1019,8 +1031,20 @@ def mergeNuclidesInto(self, compToMergeWith):
)

def iterComponents(self, typeSpec=None, exact=False):
if self.hasFlags(typeSpec, exact):
yield self
"""
Yield the component or its direct child components.

Notes
-----
This could probably be made recursive to allow arbitrary depth, but it's more complex
"""
if self._children:
for child in self._children:
if self.hasFlags(typeSpec, exact):
yield child
else:
if self.hasFlags(typeSpec, exact):
yield self
Comment on lines +1045 to +1047
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This removes the ability for component.iterComponents to yield itself. I think

if self.hasFlags(typeSpec, exact):
    yield self
# If no children, this loop isn't executed
for child in self._children:
    for c in child.iterComponents(typeSpec, exact):
        yield c

might pick up the arbitrary depth and also return the child component based on the child's flags not the flags of the parent

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

very nice, that does seem better. I'll try your suggestion.

Copy link
Member Author

@ntouran ntouran Jun 14, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok so I tried With that suggested iterComponents change, and my test_trisoCase.py unit test fails with max recursion errors.

I tried this kind of thing for a while and just couldn't quite get it working. I'm sure it can be done... it's just hard.

    [c.getNumberDensity(nuc) for nuc in nucNames]
  File "/home/nick/repos/armi/armi/reactor/composites.py", line 1256, in <listcomp>
    [c.getNumberDensity(nuc) for nuc in nucNames]
  File "/home/nick/repos/armi/armi/reactor/composites.py", line 1233, in getNumberDensity
    return self.getNuclideNumberDensities([nucName])[0]
  File "/home/nick/repos/armi/armi/reactor/components/component.py", line 605, in getNuclideNumberDensities
    composites.Composite.getNuclideNumberDensities(self, nucNames),
  File "/home/nick/repos/armi/armi/reactor/composites.py", line 1255, in getNuclideNumberDensities
    [
  File "/home/nick/repos/armi/armi/reactor/composites.py", line 1256, in <listcomp>
    [c.getNumberDensity(nuc) for nuc in nucNames]
  File "/home/nick/repos/armi/armi/reactor/composites.py", line 1256, in <listcomp>
    [c.getNumberDensity(nuc) for nuc in nucNames]
  File "/home/nick/repos/armi/armi/reactor/composites.py", line 1233, in getNumberDensity
    return self.getNuclideNumberDensities([nucName])[0]
  File "/home/nick/repos/armi/armi/reactor/components/component.py", line 605, in getNuclideNumberDensities
    composites.Composite.getNuclideNumberDensities(self, nucNames),
  File "/home/nick/repos/armi/armi/reactor/composites.py", line 1244, in getNuclideNumberDensities
    [
  File "/home/nick/repos/armi/armi/reactor/composites.py", line 1245, in <listcomp>
    c.getVolume() / (c.parent.getSymmetryFactor() if c.parent else 1.0)
  File "/home/nick/repos/armi/armi/reactor/components/component.py", line 430, in getVolume
    if self.p.volume is None:
  File "/home/nick/repos/armi/armi/reactor/parameters/parameterDefinitions.py", line 293, in __get__
    return self._getter(obj)
RecursionError: maximum recursion depth exceeded

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ooooo that's rough.

Looks like the problem could stem from Composite.getNumberDensity relying on Composite.getNuclideNumberDensities?

Comment on lines +1041 to +1047
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ntouran we've found that the mass of the matrix material is not accounted for when doing getMass calls. I think due to this section here

We modified the triso mass test blueprint file such that each element only appears in one Component. Then, block.getMass("C12") (in theory) picks up mass from the triso compact component and should be non-zero. Unfortunately this is not the case 😬

# A full blueprint file that contains assems made with ComponentGroups
blocks:
    fuel2: &block_fuel2
        triso compact: 
            shape: Circle
            blends: 
              - triso particle
            blendFracs: 
              - 0.70
            material: Graphite
            Tinput: 600.0
            Thot: 600.0
            mult: 1
            id: 0.0
            od: 1.2446
        duct: 
            shape: Hexagon
            material: Cu
            Tinput: 600.0
            Thot: 600.0
            ip: 9.0
            mult: 1.0
            op: 10.0
        coolant: 
            shape: DerivedShape
            material: NaCl
            Tinput: 25.0
            Thot: 600.0

components:
    # triso dims from Table 4 in TRISO benchmark
    kernel:
        shape: Sphere
        material: Lead
        Tinput: 600.0
        Thot: 600.0
        id: 0.0
        mult: 1.0
        od: 0.010625
    buffer:
        shape: Sphere
        material: Lead
        Tinput: 600.0
        Thot: 600.0
        id: kernel.od
        mult: 1.0
        od: 0.015625
    ipyc:
        shape: Sphere
        material: Lead
        Tinput: 600.0
        Thot: 600.0
        id: buffer.od
        mult: 1.0
        od: 0.017375
    sic:
        shape: Sphere
        material: Lead
        Tinput: 600.0
        Thot: 600.0
        id: ipyc.od
        mult: 1.0
        od: 0.019125
    opyc:
        shape: Sphere
        material: Lead
        Tinput: 600.0
        Thot: 600.0
        id: sic.od
        mult: 1.0
        od: 0.021125

component groups:
    triso particle:
      kernel:
        mult: 1
      buffer:
        mult: 1
      ipyc:
        mult: 1
      sic:
        mult: 1
      opyc:
        mult: 1


assemblies:
    fuel c: &assembly_c
        specifier: OC
        blocks: [*block_fuel2]
        height: [1.0]
        axial mesh points: [1]
        xs types: [A]
systems:
    core:
        grid name: core
        origin:
            x: 0.0
            y: 0.0
            z: 0.0
grids:
    core:
        geom: hex
        symmetry: full
        lattice map: |
          OC


nuclide flags:
    B: &nuc_flags {burn: false, xs: true}
    C: *nuc_flags
    PB: *nuc_flags
    CU63: *nuc_flags
    CU65: *nuc_flags
    CL35: *nuc_flags
    CL37: *nuc_flags
    CL37: *nuc_flags
    NA23: *nuc_flags

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @drewj-usnctech I will bring in your modified blueprints file and update the unit test so it fails as you've found and then try to fix.


def backUp(self):
"""
Expand Down
12 changes: 11 additions & 1 deletion armi/reactor/composites.py
Original file line number Diff line number Diff line change
Expand Up @@ -1237,7 +1237,13 @@ def getNumberDensity(self, nucName):
return self.getNuclideNumberDensities([nucName])[0]

def getNuclideNumberDensities(self, nucNames):
"""Return a list of number densities in atoms/barn-cm for the nuc names requested."""
"""
Return a list of number densities in atoms/barn-cm for the nuc names requested.

See Also
--------
getNumberDensities: Gets all number densities as a dict rather than a list
"""
volumes = numpy.array(
[
c.getVolume() / (c.parent.getSymmetryFactor() if c.parent else 1.0)
Expand Down Expand Up @@ -1284,6 +1290,10 @@ def getNumberDensities(self, expandFissionProducts=False):
-------
numberDensities : dict
nucName keys, number density values (atoms/bn-cm)

See Also
--------
getNuclideNumberDensities: Gets a list of number densities based on a list of nuclide names
"""
numberDensities = self._getNdensHelper()
if expandFissionProducts:
Expand Down
Loading