Skip to content

Commit

Permalink
Fix numerical diffusion in uniformMesh converter (#992)
Browse files Browse the repository at this point in the history
Mapping anything from a non-uniform assembly to a uniform assembly and
then back will cause numerical diffusion of that value. This could be a
block parameter (like mgFlux, detailedDpa, etc.) or number densities.

This commit avoids numerical diffusion by restricting the parameter mapping 
between the non-uniform and uniform assembly to only the parameters that 
are necessary inputs or designated outputs of interest for a given physics 
solver using the uniform mesh converter.

Co-authored-by: jstilley <jstilley@terrapower.com>
  • Loading branch information
mgjarrett and john-science committed Dec 14, 2022
1 parent dbc2b0c commit 0d039f2
Show file tree
Hide file tree
Showing 7 changed files with 435 additions and 126 deletions.
6 changes: 3 additions & 3 deletions armi/physics/neutronics/globalFlux/globalFluxInterface.py
Expand Up @@ -39,7 +39,6 @@
RX_ABS_MICRO_LABELS = ["nGamma", "fission", "nalph", "np", "nd", "nt"]
RX_PARAM_NAMES = ["rateCap", "rateFis", "rateProdN2n", "rateProdFis", "rateAbs"]


# pylint: disable=too-many-public-methods
class GlobalFluxInterface(interfaces.Interface):
"""
Expand Down Expand Up @@ -286,7 +285,7 @@ def __init__(self, label: Optional[str] = None):
self.real = True
self.adjoint = False
self.neutrons = True
self.photons = None
self.photons = False
self.boundaryConditions = {}
self.epsFissionSourceAvg = None
self.epsFissionSourcePoint = None
Expand Down Expand Up @@ -422,7 +421,8 @@ def _performGeometryTransformations(self, makePlots=False):
converter = self.geomConverters.get("axial")
if not converter:
if self.options.detailedAxialExpansion or self.options.hasNonUniformAssems:
converter = uniformMesh.NeutronicsUniformMeshConverter(
converterCls = uniformMesh.converterFactory(self.options)
converter = converterCls(
cs=self.options.cs,
calcReactionRates=self.options.calcReactionRatesOnMeshConversion,
)
Expand Down
Expand Up @@ -77,6 +77,18 @@ def getExecuterCls(self):
return MockGlobalFluxExecuter


class MockGlobalFluxWithExecutersNonUniform(MockGlobalFluxWithExecuters):
def getExecuterOptions(self, label=None):
"""
Return modified executerOptions
"""
opts = globalFluxInterface.GlobalFluxInterfaceUsingExecuters.getExecuterOptions(
self, label=label
)
opts.hasNonUniformAssems = True # to increase test coverage
return opts


class MockGlobalFluxExecuter(globalFluxInterface.GlobalFluxExecuter):
"""Tests for code that uses Executers, which rely on OutputReaders to update state."""

Expand Down Expand Up @@ -180,6 +192,33 @@ def test_getExecuterCls(self):
self.assertEqual(class0, globalFluxInterface.GlobalFluxExecuter)


class TestGlobalFluxInterfaceWithExecutersNonUniform(unittest.TestCase):
"""Tests for global flux execution with non-uniform assemblies."""

@classmethod
def setUpClass(cls):
cs = settings.Settings()
_o, cls.r = test_reactors.loadTestReactor()
cls.r.core.p.keff = 1.0
cls.gfi = MockGlobalFluxWithExecutersNonUniform(cls.r, cs)

def test_executerInteractionNonUniformAssems(self):
gfi, r = self.gfi, self.r
gfi.interactBOC()
gfi.interactEveryNode(0, 0)
r.p.timeNode += 1
gfi.interactEveryNode(0, 1)
gfi.interactEOC()
self.assertAlmostEqual(r.core.p.rxSwing, (1.02 - 1.01) / 1.01 * 1e5)

def test_calculateKeff(self):
self.assertEqual(self.gfi.calculateKeff(), 1.05) # set in mock

def test_getExecuterCls(self):
class0 = globalFluxInterface.GlobalFluxInterfaceUsingExecuters.getExecuterCls()
self.assertEqual(class0, globalFluxInterface.GlobalFluxExecuter)


class TestGlobalFluxResultMapper(unittest.TestCase):
"""
Test that global flux result mappings run.
Expand Down
44 changes: 36 additions & 8 deletions armi/physics/neutronics/parameters.py
Expand Up @@ -99,6 +99,7 @@ def mgFlux(self, value):
categories=[
parameters.Category.fluxQuantities,
parameters.Category.multiGroupQuantities,
parameters.Category.gamma,
],
default=None,
)
Expand Down Expand Up @@ -129,7 +130,10 @@ def mgFlux(self, value):
description="multigroup gamma source",
location=ParamLocation.AVERAGE,
saveToDB=True,
categories=[parameters.Category.multiGroupQuantities],
categories=[
parameters.Category.multiGroupQuantities,
parameters.Category.gamma,
],
default=None,
)

Expand All @@ -139,7 +143,8 @@ def mgFlux(self, value):
description="gamma source",
location=ParamLocation.AVERAGE,
saveToDB=True,
default=None,
categories=[parameters.Category.gamma],
default=0.0,
)

pb.defParam(
Expand Down Expand Up @@ -183,7 +188,7 @@ def mgFlux(self, value):
"pinMgFluxesGamma",
units="g/s/cm$^2$",
description="should be a blank 3-D array, but re-defined later (ng x nPins x nAxialSegments)",
categories=[parameters.Category.pinQuantities],
categories=[parameters.Category.pinQuantities, parameters.Category.gamma],
saveToDB=False,
default=None,
)
Expand Down Expand Up @@ -289,13 +294,14 @@ def linPowByPinNeutron(self, value):
else:
self._p_linPowByPinNeutron = numpy.array(value)

# gamma category because linPowByPin is only split by neutron/gamma when gamma is activated
pb.defParam(
"linPowByPinNeutron",
setter=linPowByPinNeutron,
units="W/cm",
description="Pin linear neutron heat rate. This is the neutron heating component of `linPowByPin`",
location=ParamLocation.CHILDREN,
categories=[parameters.Category.pinQuantities],
categories=[parameters.Category.pinQuantities, parameters.Category.gamma],
default=None,
)

Expand All @@ -311,7 +317,7 @@ def linPowByPinGamma(self, value):
units="W/cm",
description="Pin linear gamma heat rate. This is the gamma heating component of `linPowByPin`",
location=ParamLocation.CHILDREN,
categories=[parameters.Category.pinQuantities],
categories=[parameters.Category.pinQuantities, parameters.Category.gamma],
default=None,
)

Expand All @@ -320,6 +326,7 @@ def linPowByPinGamma(self, value):
units="#/s",
description='List of reaction rates in specified by setting "reactionsToDB"',
location=ParamLocation.VOLUME_INTEGRATED,
categories=[parameters.Category.fluxQuantities],
default=None,
)

Expand Down Expand Up @@ -546,17 +553,25 @@ def linPowByPinGamma(self, value):
"pdensGamma",
units="W/cm^3",
description="Average volumetric gamma power density",
categories=[parameters.Category.gamma],
)

# gamma category because pdens is only split by neutron/gamma when gamma is activated
pb.defParam(
"pdensNeutron",
units="W/cm^3",
description="Average volumetric neutron power density",
categories=[parameters.Category.gamma],
)

pb.defParam("ppdens", units="W/cm^3", description="Peak power density")

pb.defParam("ppdensGamma", units="W/cm^3", description="Peak gamma density")
pb.defParam(
"ppdensGamma",
units="W/cm^3",
description="Peak gamma density",
categories=[parameters.Category.gamma],
)

# rx rate params that are derived during mesh conversion.
# We'd like all things that can be derived from flux and XS to be
Expand Down Expand Up @@ -603,15 +618,27 @@ def linPowByPinGamma(self, value):
"powerGenerated",
units=" W",
description="Generated power. Different than b.p.power only when gamma transport is activated.",
categories=[parameters.Category.gamma],
)

pb.defParam("power", units="W", description="Total power")

pb.defParam("powerDecay", units="W", description="Total decay power")

pb.defParam("powerGamma", units="W", description="Total gamma power")
pb.defParam(
"powerGamma",
units="W",
description="Total gamma power",
categories=[parameters.Category.gamma],
)

pb.defParam("powerNeutron", units="W", description="Total neutron power")
# gamma category because power is only split by neutron/gamma when gamma is activated
pb.defParam(
"powerNeutron",
units="W",
description="Total neutron power",
categories=[parameters.Category.gamma],
)

with pDefs.createBuilder(default=0.0) as pb:
pb.defParam(
Expand Down Expand Up @@ -663,6 +690,7 @@ def linPowByPinGamma(self, value):
units="W/cm^3",
description="Volume-averaged generated power density. Different than b.p.pdens only when gamma transport is activated.",
location=ParamLocation.AVERAGE,
categories=[parameters.Category.gamma],
)

return pDefs
Expand Down
130 changes: 128 additions & 2 deletions armi/reactor/converters/tests/test_uniformMesh.py
Expand Up @@ -28,6 +28,28 @@
from armi.tests import TEST_ROOT, ISOAA_PATH


class DummyFluxOptions:
def __init__(self):
self.photons = False


class TestConverterFactory(unittest.TestCase):
def setUp(self):
self.o, self.r = test_reactors.loadTestReactor(
inputFilePath=os.path.join(TEST_ROOT, "detailedAxialExpansion"),
)
self.dummyOptions = DummyFluxOptions()

def test_converterFactory(self):
self.dummyOptions.photons = False
neutronConverter = uniformMesh.converterFactory(self.dummyOptions)
self.assertTrue(neutronConverter, uniformMesh.NeutronicsUniformMeshConverter)

self.dummyOptions.photons = True
gammaConverter = uniformMesh.converterFactory(self.dummyOptions)
self.assertTrue(gammaConverter, uniformMesh.GammaUniformMeshConverter)


class TestAssemblyUniformMesh(unittest.TestCase):
"""
Tests individual operations of the uniform mesh converter
Expand All @@ -48,7 +70,10 @@ def test_makeAssemWithUniformMesh(self):

self.converter._computeAverageAxialMesh()
newAssem = self.converter.makeAssemWithUniformMesh(
sourceAssem, self.converter._uniformMesh
sourceAssem,
self.converter._uniformMesh,
blockParamNames=["power"],
mapNumberDensities=True,
)

prevB = None
Expand Down Expand Up @@ -89,7 +114,7 @@ def test_makeAssemWithUniformMeshSubmesh(self):

self.r.core.updateAxialMesh()
newAssem = self.converter.makeAssemWithUniformMesh(
sourceAssem, self.r.core.p.axialMesh[1:]
sourceAssem, self.r.core.p.axialMesh[1:], blockParamNames=["power"]
)

self.assertNotEqual(len(newAssem), len(sourceAssem))
Expand All @@ -114,6 +139,7 @@ def test_makeAssemUniformMeshParamMappingSameMesh(self):
sourceAssem,
sourceAssem.getAxialMesh(),
blockParamNames=["flux", "power", "mgFlux"],
mapNumberDensities=True,
)
for b, origB in zip(newAssem, sourceAssem):
self.assertEqual(b.p.flux, 1.0)
Expand Down Expand Up @@ -329,6 +355,106 @@ def test_applyStateToOriginal(self):
)


class TestGammaUniformMesh(unittest.TestCase):
"""
Tests gamma uniform mesh converter
Loads reactor once per test
"""

@classmethod
def setUpClass(cls):
# random seed to support random mesh in unit tests below
random.seed(987324987234)

def setUp(self):
self.o, self.r = test_reactors.loadTestReactor(
TEST_ROOT, customSettings={"xsKernel": "MC2v2"}
)
self.r.core.lib = isotxs.readBinary(ISOAA_PATH)
self.r.core.p.keff = 1.0
self.converter = uniformMesh.GammaUniformMeshConverter(
cs=self.o.cs, calcReactionRates=False
)

def test_convertNumberDensities(self):
refMass = self.r.core.getMass("U235")
applyNonUniformHeightDistribution(
self.r
) # this changes the mass of everything in the core
perturbedCoreMass = self.r.core.getMass("U235")
self.assertNotEqual(refMass, perturbedCoreMass)
self.converter.convert(self.r)

uniformReactor = self.converter.convReactor
uniformMass = uniformReactor.core.getMass("U235")

self.assertAlmostEqual(
perturbedCoreMass, uniformMass
) # conversion conserved mass
self.assertAlmostEqual(
self.r.core.getMass("U235"), perturbedCoreMass
) # conversion didn't change source reactor mass

def test_applyStateToOriginal(self):
applyNonUniformHeightDistribution(self.r) # note: this perturbs the ref. mass

# set original parameters on pre-mapped core with non-uniform assemblies
for b in self.r.core.getBlocks():
b.p.mgFlux = range(33)
b.p.adjMgFlux = range(33)
b.p.fastFlux = 2.0
b.p.flux = 5.0
b.p.power = 5.0

# set original parameters on pre-mapped core with non-uniform assemblies
self.converter.convert(self.r)
for b in self.converter.convReactor.core.getBlocks():
b.p.powerGamma = 0.5
b.p.powerNeutron = 0.5

# check integral and density params
assemblyPowers = [
a.calcTotalParam("power") for a in self.converter.convReactor.core
]
assemblyGammaPowers = [
a.calcTotalParam("powerGamma") for a in self.converter.convReactor.core
]
totalPower = self.converter.convReactor.core.calcTotalParam(
"power", generationNum=2
)
totalPowerGamma = self.converter.convReactor.core.calcTotalParam(
"powerGamma", generationNum=2
)

self.converter.applyStateToOriginal()

for b in self.r.core.getBlocks():
# equal to original value because these were never mapped
self.assertEqual(b.p.fastFlux, 2.0)
self.assertEqual(b.p.flux, 5.0)
self.assertEqual(b.p.fastFlux, 2.0)
self.assertEqual(b.p.power, 5.0)

# not equal because blocks are different size
self.assertNotEqual(b.p.powerGamma, 0.5)
self.assertNotEqual(b.p.powerNeutron, 0.5)

# equal because these are mapped
for expectedPower, expectedGammaPower, a in zip(
assemblyPowers, assemblyGammaPowers, self.r.core
):
self.assertAlmostEqual(a.calcTotalParam("power"), expectedPower)
self.assertAlmostEqual(a.calcTotalParam("powerGamma"), expectedGammaPower)

self.assertAlmostEqual(
self.r.core.calcTotalParam("powerGamma", generationNum=2), totalPowerGamma
)
self.assertAlmostEqual(
self.r.core.calcTotalParam("power", generationNum=2), totalPower
)


class TestParamConversion(unittest.TestCase):
def setUp(self):
"""
Expand Down

0 comments on commit 0d039f2

Please sign in to comment.