Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/main' into drewj/improve-assem-a…
Browse files Browse the repository at this point in the history 
…xial-linkage

* origin/main:
  New plugin hook for before reactor construction (#1945)
  Removing DoseResultsMapper (#1952)
  Adding logic so HexAssemblies are hexagonal (#1935)
  Adding Block.getInputHeight (#1927)
  Removing Assembly.doubleResolution (#1951)
  Adding DeprecationWarning for HistoryTrackerInterface (#1950)
  Removing tabulate as a dependency (#1948)
  • Loading branch information
@drewj-tp
drewj-tp committed Oct 16, 2024
2 parents 83ee423 + 01b0ade commit d925b98
Show file tree
Hide file tree
Showing 14 changed files with 191 additions and 540 deletions.
3 changes: 3 additions & 0 deletions armi/bookkeeping/historyTracker.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,6 +139,9 @@ def __init__(self, r, cs):
self.xsHistory = {}
self._preloadedBlockHistory = None

msg = "The HistoryTrackerInterface is deprecated, and will be removed."
runLog.warning(msg)

def interactBOL(self):
self.addDetailAssembliesBOL()

Expand Down
336 changes: 0 additions & 336 deletions armi/physics/neutronics/globalFlux/globalFluxInterface.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,6 @@
from typing import Dict, Optional

import numpy as np
import scipy.integrate

from armi import interfaces
from armi import runLog
Expand Down Expand Up @@ -894,341 +893,6 @@ def _updateAssemblyLevelParams(self):
a.p.kInf = totalSrc / totalAbs # assembly average k-inf.


class DoseResultsMapper(GlobalFluxResultMapper):
"""
Updates fluence and dpa when time shifts.
Often called after a depletion step. It is invoked using :py:meth:`apply() <.DoseResultsMapper.apply>`.
Parameters
----------
depletionSeconds: float, required
Length of depletion step in units of seconds
options: GlobalFluxOptions, required
Object describing options used by the global flux solver. A few attributes from
this object are used to run the methods in DoseResultsMapper. An example
attribute is aclpDoseLimit.
Notes
-----
We attempted to make this a set of stateless functions but the requirement of various
options made it more of a data passing task than we liked. So it's just a lightweight
and ephemeral results mapper.
"""

def __init__(self, depletionSeconds, options):
self.success = False
self.options = options
self.cs = self.options.cs
self.r = None
self.depletionSeconds = depletionSeconds

def apply(self, reactor, blockList=None):
"""
Invokes :py:meth:`updateFluenceAndDpa() <.DoseResultsMapper.updateFluenceAndDpa>`
for a provided Reactor object.
Parameters
----------
reactor: Reactor, required
ARMI Reactor object
blockList: list, optional
List of ARMI blocks to be processed by the class. If no blocks are provided, then
blocks returned by :py:meth:`getBlocks() <.reactors.Core.getBlocks>` are used.
Returns
-------
None
"""
runLog.extra("Updating fluence and dpa on reactor based on depletion step.")
self.r = reactor
self.updateFluenceAndDpa(self.depletionSeconds, blockList=blockList)

def updateFluenceAndDpa(self, stepTimeInSeconds, blockList=None):
r"""
Updates the fast fluence and the DPA of the blocklist.
The dpa rate from the previous timestep is used to compute the dpa here.
There are several items of interest computed here, including:
* detailedDpa: The average DPA of a block
* detailedDpaPeak: The peak dpa of a block, considering axial and radial peaking
The peaking is based either on a user-provided peaking factor (computed in a
pin reconstructed rotation study) or the nodal flux peaking factors
* dpaPeakFromFluence: fast fluence * fluence conversion factor (old and inaccurate).
Used to be dpaPeak
Parameters
----------
stepTimeInSeconds : float
Time in seconds that the cycle ran at the current mgFlux
blockList : list, optional
blocks to be updated. Defaults to all blocks in core
See Also
--------
updateDpaRate : updates the DPA rate used here to compute actual dpa
"""
blockList = blockList or self.r.core.getBlocks()

if not blockList[0].p.fluxPeak:
runLog.warning(
"no fluxPeak parameter on this model. Peak DPA will be equal to average DPA. "
"Perhaps you are not running a nodal approximation."
)

for b in blockList:
b.p.residence += stepTimeInSeconds / units.SECONDS_PER_DAY
b.p.fluence += b.p.flux * stepTimeInSeconds
b.p.fastFluence += b.p.flux * stepTimeInSeconds * b.p.fastFluxFr
b.p.fastFluencePeak += b.p.fluxPeak * stepTimeInSeconds * b.p.fastFluxFr

# update detailed DPA based on dpa rate computed at LAST timestep.
# new incremental DPA increase for duct distortion interface (and eq)
b.p.newDPA = b.p.detailedDpaRate * stepTimeInSeconds
b.p.newDPAPeak = b.p.detailedDpaPeakRate * stepTimeInSeconds

# use = here instead of += because we need the param system to notice the change for syncronization.
b.p.detailedDpa = b.p.detailedDpa + b.p.newDPA
b.p.detailedDpaPeak = b.p.detailedDpaPeak + b.p.newDPAPeak
b.p.detailedDpaThisCycle = b.p.detailedDpaThisCycle + b.p.newDPA

# increment point dpas
# this is specific to hex geometry, but they are general neutronics block parameters
# if it is a non-hex block, this should be a no-op
if b.p.pointsCornerDpaRate is not None:
if b.p.pointsCornerDpa is None:
b.p.pointsCornerDpa = np.zeros((6,))
b.p.pointsCornerDpa = (
b.p.pointsCornerDpa + b.p.pointsCornerDpaRate * stepTimeInSeconds
)
if b.p.pointsEdgeDpaRate is not None:
if b.p.pointsEdgeDpa is None:
b.p.pointsEdgeDpa = np.zeros((6,))
b.p.pointsEdgeDpa = (
b.p.pointsEdgeDpa + b.p.pointsEdgeDpaRate * stepTimeInSeconds
)

if self.options.dpaPerFluence:
# do the less rigorous fluence -> DPA conversion if the user gave a factor.
b.p.dpaPeakFromFluence = (
b.p.fastFluencePeak * self.options.dpaPerFluence
)

# Set burnup peaking
# b.p.percentBu/buRatePeak is expected to have been updated elsewhere (depletion)
# (this should run AFTER burnup has been updated)
# try to find the peak rate first
peakRate = None
if b.p.buRatePeak:
# best case scenario, we have peak burnup rate
peakRate = b.p.buRatePeak
elif b.p.buRate:
# use whatever peaking factor we can find if just have rate
peakRate = b.p.buRate * self.getBurnupPeakingFactor(b)

# If peak rate found, use to calc peak burnup; otherwise scale burnup
if peakRate:
# peakRate is in per day
peakRatePerSecond = peakRate / units.SECONDS_PER_DAY
b.p.percentBuPeak = (
b.p.percentBuPeak + peakRatePerSecond * stepTimeInSeconds
)
else:
# No rate, make bad assumption.... assumes peaking is same at each position through
# shuffling/irradiation history...
runLog.warning(
"Scaling burnup by current peaking factor... This assumes peaking "
"factor was constant through shuffling/irradiation history.",
single=True,
)
b.p.percentBuPeak = b.p.percentBu * self.getBurnupPeakingFactor(b)

for a in self.r.core.getAssemblies():
a.p.daysSinceLastMove += stepTimeInSeconds / units.SECONDS_PER_DAY

self.updateMaxDpaParams()
self.updateCycleDoseParams()
self.updateLoadpadDose()

def updateCycleDoseParams(self):
"""Updates reactor params based on the amount of dose (detailedDpa) accrued this cycle.
Params updated include:
* maxDetailedDpaThisCycle
* dpaFullWidthHalfMax
* elevationOfACLP3Cycles
* elevationOfACLP7Cycles
These parameters are left as zeroes at BOC because no dose has been accumulated yet.
"""
cycle = self.r.p.cycle
timeNode = self.r.p.timeNode

if timeNode <= 0:
return

daysIntoCycle = sum(self.r.o.stepLengths[cycle][:timeNode])
cycleLength = self.r.p.cycleLength

maxDetailedDpaThisCycle = 0.0
peakDoseAssem = None
for a in self.r.core:
if a.getMaxParam("detailedDpaThisCycle") > maxDetailedDpaThisCycle:
maxDetailedDpaThisCycle = a.getMaxParam("detailedDpaThisCycle")
peakDoseAssem = a
self.r.core.p.maxDetailedDpaThisCycle = maxDetailedDpaThisCycle

if peakDoseAssem is None:
return

doseHalfMaxHeights = peakDoseAssem.getElevationsMatchingParamValue(
"detailedDpaThisCycle", maxDetailedDpaThisCycle / 2.0
)
if len(doseHalfMaxHeights) != 2:
runLog.warning(
"Something strange with detailedDpaThisCycle shape in {}, "
"non-2 number of values matching {}".format(
peakDoseAssem, maxDetailedDpaThisCycle / 2.0
)
)
else:
self.r.core.p.dpaFullWidthHalfMax = (
doseHalfMaxHeights[1] - doseHalfMaxHeights[0]
)

aclpDoseLimit = self.options.aclpDoseLimit
aclpDoseLimit3 = aclpDoseLimit / 3.0 * (daysIntoCycle / cycleLength)
aclpLocations3 = peakDoseAssem.getElevationsMatchingParamValue(
"detailedDpaThisCycle", aclpDoseLimit3
)
if len(aclpLocations3) != 2:
runLog.warning(
"Something strange with detailedDpaThisCycle shape in {}"
", non-2 number of values matching {}".format(
peakDoseAssem, aclpDoseLimit / 3.0
)
)
else:
self.r.core.p.elevationOfACLP3Cycles = aclpLocations3[1]

aclpDoseLimit7 = aclpDoseLimit / 7.0 * (daysIntoCycle / cycleLength)
aclpLocations7 = peakDoseAssem.getElevationsMatchingParamValue(
"detailedDpaThisCycle", aclpDoseLimit7
)
if len(aclpLocations7) != 2:
runLog.warning(
"Something strange with detailedDpaThisCycle shape in {}, "
"non-2 number of values matching {}".format(
peakDoseAssem, aclpDoseLimit / 7.0
)
)
else:
self.r.core.p.elevationOfACLP7Cycles = aclpLocations7[1]

def updateLoadpadDose(self):
"""
Summarize the dose in DPA of the above-core load pad.
This sets the following reactor params:
* loadPadDpaPeak : the peak dpa in any load pad
* loadPadDpaAvg : the highest average dpa in any load pad
.. warning:: This only makes sense for cores with load
pads on their assemblies.
See Also
--------
_calcLoadPadDose : computes the load pad dose
"""
peakPeak, peakAvg = self._calcLoadPadDose()
if peakPeak is None:
return
self.r.core.p.loadPadDpaAvg = peakAvg[0]
self.r.core.p.loadPadDpaPeak = peakPeak[0]
str_ = [
"Above-core load pad (ACLP) summary. It starts at {0} cm above the "
"bottom of the grid plate".format(self.options.loadPadElevation)
]
str_.append(
"The peak ACLP dose is {0} DPA in {1}".format(peakPeak[0], peakPeak[1])
)
str_.append(
"The max avg. ACLP dose is {0} DPA in {1}".format(peakAvg[0], peakAvg[1])
)
runLog.info("\n".join(str_))

def _calcLoadPadDose(self):
r"""
Determines some dose information on the load pads if they exist.
Expects a few settings to be present in cs
loadPadElevation : float
The bottom of the load pad's elevation in cm above the bottom of
the (upper) grid plate.
loadPadLength : float
The axial length of the load pad to average over
This builds axial splines over the assemblies and then integrates them over the load pad.
The assumptions are that detailedDpa is the average, defined in the center
and detailedDpaPeak is the peak, also defined in the center of blocks.
Parameters
----------
core : armi.reactor.reactors.Core
cs : armi.settings.caseSettings.Settings
Returns
-------
peakPeak : tuple
A (peakValue, peakAssem) tuple
peakAvg : tuple
a (avgValue, avgAssem) tuple
See Also
--------
writeLoadPadDoseSummary : prints out the dose
Assembly.getParamValuesAtZ : gets the parameters at any arbitrary z point
"""
loadPadBottom = self.options.loadPadElevation
loadPadLength = self.options.loadPadLength
if not loadPadBottom or not loadPadLength:
# no load pad dose requested
return None, None

peakPeak = (0.0, None)
peakAvg = (0.0, None)
loadPadTop = loadPadBottom + loadPadLength

zrange = np.linspace(loadPadBottom, loadPadTop, 100)
for a in self.r.core.getAssemblies(Flags.FUEL):
# scan over the load pad to find the peak dpa
# no caching.
peakDose = max(
a.getParamValuesAtZ("detailedDpaPeak", zrange, fillValue="extrapolate")
)
# restrict to fuel because control assemblies go nuts in dpa.
integrand = a.getParamOfZFunction("detailedDpa", fillValue="extrapolate")
returns = scipy.integrate.quad(integrand, loadPadBottom, loadPadTop)
avgDose = (
float(returns[0]) / loadPadLength
) # convert to float in case it's an ndarray

# track max doses
if peakDose > peakPeak[0]:
peakPeak = (peakDose, a)
if avgDose > peakAvg[0]:
peakAvg = (avgDose, a)

return peakPeak, peakAvg


def computeDpaRate(mgFlux, dpaXs):
r"""
Compute the DPA rate incurred by exposure of a certain flux spectrum.
Expand Down
Loading

0 comments on commit d925b98

Please sign in to comment.