Repository: terrapower/armi

Issue:
Support specifying total reactor power through power density
For some analysis cases, it is easier to investigate concepts that share a power density but may differ in total size or fuel mass. Therefore the total power would differ. Meaning if one were to create a [parameter sweep](https://terrapower.github.io/armi/tutorials/param_sweep.html) that tried to keep power *density* constant, there are some hoops to jump through.
1. Settings modifier to modify power : pretty easy to implement once you know the power.
2. Other input modifiers to perturb your reactor

But, input modifiers work on the blueprints, so the entire reactor isn't created to be able to do something like `power = powerDensity * reactor.getVolume()` (for a W/cc power density) or `powerDensity * reactor.getHMMass()` (for W/gHM power density)

This ticket proposes allowing the total power to be set via power density specified in the settings file. Or maybe two power density settings, for volumetric and gravimetric power density)

Some settings pairs could be
- `powerPerVolume` and `powerPerHM` (short but no units explicitly stated)
- `volumetricPower` and `gravimetricPower` (no units, a little archaic naming?)
- `powerPerCC` and `powerPerGHM` (units are there, maybe confusing `GHM` for grams HM or giga-HM)

These settings should be mutually exclusive with each other and with `power`. But I think after some initial conversion, the total power should be the metric used by ARMI like it is already.

Discussion:

Patch:
diff --git a/armi/operators/operator.py b/armi/operators/operator.py
index d8c7fc5c..6196c215 100644
--- a/armi/operators/operator.py
+++ b/armi/operators/operator.py
@@ -355,10 +355,13 @@ class Operator:
         if halt:
             return False
 
+        # read total core power from settings (power or powerDensity)
+        basicPower = self.cs["power"] or (
+            self.cs["powerDensity"] * self.r.core.getHMMass()
+        )
+
         for timeNode in range(startingNode, int(self.burnSteps[cycle])):
-            self.r.core.p.power = (
-                self.powerFractions[cycle][timeNode] * self.cs["power"]
-            )
+            self.r.core.p.power = self.powerFractions[cycle][timeNode] * basicPower
             self.r.p.capacityFactor = (
                 self.r.p.availabilityFactor * self.powerFractions[cycle][timeNode]
             )
@@ -373,7 +376,7 @@ class Operator:
             else:
                 powFrac = self.powerFractions[cycle][timeNode - 1]
 
-            self.r.core.p.power = powFrac * self.cs["power"]
+            self.r.core.p.power = powFrac * basicPower
             self._timeNodeLoop(cycle, timeNode)
 
         self.interactAllEOC(self.r.p.cycle)
diff --git a/armi/operators/settingsValidation.py b/armi/operators/settingsValidation.py
index bd0fa4f0..ebab2cbc 100644
--- a/armi/operators/settingsValidation.py
+++ b/armi/operators/settingsValidation.py
@@ -20,25 +20,25 @@ These then pop up during initialization of a run, either on the command line or
 dialogues in the GUI. They say things like: "Your ___ setting has the value ___, which
 is impossible. Would you like to switch to ___?"
 """
-import re
+import itertools
 import os
+import re
 import shutil
-import itertools
 
 from armi import context
 from armi import getPluginManagerOrFail
 from armi import runLog
-from armi.utils import pathTools
-from armi.utils.mathematics import expandRepeatedFloats
+from armi.physics import neutronics
 from armi.reactor import geometry
 from armi.reactor import systemLayoutInput
-from armi.physics import neutronics
-from armi.utils import directoryChangers
 from armi.settings.settingsIO import (
     prompt,
     RunLogPromptCancel,
     RunLogPromptUnresolvable,
 )
+from armi.utils import directoryChangers
+from armi.utils import pathTools
+from armi.utils.mathematics import expandRepeatedFloats
 
 
 class Query:
@@ -468,8 +468,15 @@ class Inspector:
         )
 
         self.addQuery(
-            lambda: not self.cs["power"],
-            "No power level set. You must always start by importing a base settings file.",
+            lambda: not self.cs["power"] and not self.cs["powerDensity"],
+            "No power or powerDensity set. You must always start by importing a base settings file.",
+            "",
+            self.NO_ACTION,
+        )
+
+        self.addQuery(
+            lambda: self.cs["power"] > 0 and self.cs["powerDensity"] > 0,
+            "The power and powerDensity are both set, please note the power will be used as the truth.",
             "",
             self.NO_ACTION,
         )
diff --git a/armi/operators/snapshots.py b/armi/operators/snapshots.py
index d4337124..f5865462 100644
--- a/armi/operators/snapshots.py
+++ b/armi/operators/snapshots.py
@@ -73,6 +73,7 @@ class OperatorSnapshots(operatorMPI.OperatorMPI):
             # need to update reactor power after the database load
             # this is normally handled in operator._cycleLoop
             self.r.core.p.power = self.cs["power"]
+            self.r.core.p.powerDensity = self.cs["powerDensity"]
 
             halt = self.interactAllBOC(self.r.p.cycle)
             if halt:
diff --git a/armi/physics/neutronics/globalFlux/globalFluxInterface.py b/armi/physics/neutronics/globalFlux/globalFluxInterface.py
index 8086fad3..5c5d2047 100644
--- a/armi/physics/neutronics/globalFlux/globalFluxInterface.py
+++ b/armi/physics/neutronics/globalFlux/globalFluxInterface.py
@@ -127,6 +127,7 @@ class GlobalFluxInterface(interfaces.Interface):
             )
             / units.WATTS_PER_MW
         )
+        self.r.core.setPowerIfNecessary()
         specifiedPower = (
             self.r.core.p.power / units.WATTS_PER_MW / self.r.core.powerMultiplier
         )
diff --git a/armi/reactor/reactorParameters.py b/armi/reactor/reactorParameters.py
index c3097974..e738535e 100644
--- a/armi/reactor/reactorParameters.py
+++ b/armi/reactor/reactorParameters.py
@@ -366,6 +366,14 @@ def defineCoreParameters():
             "nuclear power generated in the core.",
         )
 
+        pb.defParam(
+            "powerDensity",
+            units=f"{units.WATTS}/{units.GRAMS}",
+            description="BOL Power density of the reactor core, in units of Watts per"
+            "grams of Heavy Metal Mass. After the BOL, the power parameter will be set, "
+            "and this will entirely overridden by that.",
+        )
+
         pb.defParam(
             "powerDecay",
             units=units.WATTS,
diff --git a/armi/reactor/reactors.py b/armi/reactor/reactors.py
index 86721e4a..4782f15a 100644
--- a/armi/reactor/reactors.py
+++ b/armi/reactor/reactors.py
@@ -395,8 +395,8 @@ class Core(composites.Composite):
         The reference assembly is defined as the center-most assembly with a FUEL flag,
         if any are present, or the center-most of any assembly otherwise.
 
-        Warnings
-        ========
+        Warning
+        -------
         The convenience of this property should be weighed against it's somewhat
         arbitrary nature for any particular client. The center-most fueled assembly is
         not particularly representative of the state of the core as a whole.
@@ -443,6 +443,20 @@ class Core(composites.Composite):
             )
         )
 
+    def setPowerFromDensity(self):
+        """Set the power from the powerDensity."""
+        self.p.power = self.p.powerDensity * self.getHMMass()
+
+    def setPowerIfNecessary(self):
+        """Set the core power, from the power density.
+
+        If the power density is set, but the power isn't, we set the calculate the
+        total heavy metal mass of the reactor, and set the total power. Which will
+        then be the real source of truth again.
+        """
+        if self.p.power == 0 and self.p.powerDensity > 0:
+            self.setPowerFromDensity()
+
     def setBlockMassParams(self):
         """Set the parameters kgHM and kgFis for each block and calculate Pu fraction."""
         for b in self.getBlocks():
@@ -733,7 +747,9 @@ class Core(composites.Composite):
         """
         Returns the number of rings in this reactor. Based on location so indexing will start at 1.
 
-        WARNING: If you loop through range(maxRing) then ring+1 is the one you want!!
+        Warning
+        -------
+        If you loop through range(maxRing) then ring+1 is the one you want!
 
         Parameters
         ----------
@@ -742,7 +758,6 @@ class Core(composites.Composite):
 
         When circular ring shuffling is activated, this changes interpretation.
         Developers plan on making this another method for the secondary interpretation.
-
         """
         if self.circularRingList and not indexBased:
             return max(self.circularRingList)
@@ -1299,8 +1314,8 @@ class Core(composites.Composite):
         """
         Gets the first assembly in the reactor.
 
-        Warnings
-        --------
+        Warning
+        -------
         This function should be used with great care. There are **very** few
         circumstances in which one wants the "first" of a given sort of assembly,
         `whichever that may happen to be`. Precisely which assembly is returned is
diff --git a/armi/settings/fwSettings/globalSettings.py b/armi/settings/fwSettings/globalSettings.py
index 3bf5cfb8..72410a04 100644
--- a/armi/settings/fwSettings/globalSettings.py
+++ b/armi/settings/fwSettings/globalSettings.py
@@ -91,6 +91,7 @@ CONF_OUTPUT_FILE_EXTENSION = "outputFileExtension"
 CONF_PHYSICS_FILES = "savePhysicsFiles"
 CONF_PLOTS = "plots"
 CONF_POWER = "power"
+CONF_POWER_DENSITY = "powerDensity"
 CONF_POWER_FRACTIONS = "powerFractions"
 CONF_PROFILE = "profile"
 CONF_REALLY_SMALL_RUN = "reallySmallRun"
@@ -627,6 +628,14 @@ def defineSettings() -> List[setting.Setting]:
             "setting the powerFractions setting.",
             schema=vol.All(vol.Coerce(float), vol.Range(min=0)),
         ),
+        setting.Setting(
+            CONF_POWER_DENSITY,
+            default=0.0,
+            label="Reactor Thermal Power Density (W/HMM)",
+            description="Thermal power of the Reactor, per gram of Heavy metal "
+            "mass. Ignore this setting if the `power` setting is non-zero.",
+            schema=vol.All(vol.Coerce(float), vol.Range(min=0)),
+        ),
         setting.Setting(
             CONF_REMOVE_PER_CYCLE, default=3, label="Move per Cycle", description="None"
         ),
diff --git a/doc/release/0.2.rst b/doc/release/0.2.rst
index e093ddae..d04ee417 100644
--- a/doc/release/0.2.rst
+++ b/doc/release/0.2.rst
@@ -12,6 +12,7 @@ What's new in ARMI
 #. Broad cleanup of ``Parameters``: filled in all empty units and descriptions, removed unused params. (`PR#1345 <https://github.com/terrapower/armi/pull/1345>`_)
 #. Removed redundant ``Material.name`` variable. (`PR#1335 <https://github.com/terrapower/armi/pull/1335>`_)
 #. Moved the ``Reactor`` assembly number from the global scope to a ``Parameter``. (`PR#1383 <https://github.com/terrapower/armi/pull/1383>`_)
+#. Added ``powerDensity`` as a high-level alternative to ``power`` to configure a Reactor. (`PR#1395 <https://github.com/terrapower/armi/pull/1395>`_)
 #. Added SHA1 hashes of XS control files to the welcome text. (`PR#1334 <https://github.com/terrapower/armi/pull/1334>`_)
 #. Add python 3.11 to ARMI's CI testing GH actions! (`PR#1341 <https://github.com/terrapower/armi/pull/1341>`_)
 #. Put back ``avgFuelTemp`` block parameter. (`PR#1362 <https://github.com/terrapower/armi/pull/1362>`_)


Test patch:
diff --git a/armi/bookkeeping/db/tests/test_databaseInterface.py b/armi/bookkeeping/db/tests/test_databaseInterface.py
index 9570a1b4..2b653b1a 100644
--- a/armi/bookkeeping/db/tests/test_databaseInterface.py
+++ b/armi/bookkeeping/db/tests/test_databaseInterface.py
@@ -129,6 +129,7 @@ class TestDatabaseWriter(unittest.TestCase):
         self.td = directoryChangers.TemporaryDirectoryChanger()
         self.td.__enter__()
         cs = settings.Settings(os.path.join(TEST_ROOT, "armiRun.yaml"))
+        cs = cs.modified(newSettings={"power": 0.0, "powerDensity": 9e4})
         self.o, cs = getSimpleDBOperator(cs)
         self.r = self.o.r
 
@@ -139,6 +140,9 @@ class TestDatabaseWriter(unittest.TestCase):
         def goodMethod(cycle, node):
             pass
 
+        # the power should start at zero
+        self.assertEqual(self.r.core.p.power, 0)
+
         self.o.interfaces.append(MockInterface(self.o.r, self.o.cs, goodMethod))
         with self.o:
             self.o.operate()
@@ -159,6 +163,9 @@ class TestDatabaseWriter(unittest.TestCase):
             self.assertIn("settings", h5["inputs"])
             self.assertIn("blueprints", h5["inputs"])
 
+        # after operating, the power will be greater than zero
+        self.assertGreater(self.r.core.p.power, 1e9)
+
     def test_metaDataEndFail(self):
         def failMethod(cycle, node):
             if cycle == 0 and node == 1:
diff --git a/armi/conftest.py b/armi/conftest.py
index 1c1e519d..7eb1a52e 100644
--- a/armi/conftest.py
+++ b/armi/conftest.py
@@ -21,7 +21,7 @@ when using the ARMI framework as a submodule in a larger project.
 
 Tests must be invoked via pytest for this to have any affect, for example::
 
-    $ pytest -n6 framework/armi
+    $ pytest -n 6 armi
 
 """
 import os
diff --git a/armi/reactor/tests/test_reactors.py b/armi/reactor/tests/test_reactors.py
index c6e2e713..43056fbd 100644
--- a/armi/reactor/tests/test_reactors.py
+++ b/armi/reactor/tests/test_reactors.py
@@ -1098,7 +1098,7 @@ class HexReactorTests(ReactorTests):
         self.assertEqual(len(list(self.r.core.zones)), 0)
 
     def test_getNuclideCategories(self):
-        # test that nuclides are categorized correctly
+        """Test that nuclides are categorized correctly."""
         self.r.core.getNuclideCategories()
         self.assertIn("coolant", self.r.core._nuclideCategories)
         self.assertIn("structure", self.r.core._nuclideCategories)
@@ -1107,6 +1107,27 @@ class HexReactorTests(ReactorTests):
         self.assertIn("FE56", self.r.core._nuclideCategories["structure"])
         self.assertIn("U235", self.r.core._nuclideCategories["fuel"])
 
+    def test_setPowerIfNecessary(self):
+        self.assertAlmostEqual(self.r.core.p.power, 0)
+        self.assertAlmostEqual(self.r.core.p.powerDensity, 0)
+
+        # to start, this method shouldn't do anything
+        self.r.core.setPowerIfNecessary()
+        self.assertAlmostEqual(self.r.core.p.power, 0)
+
+        # take the powerDensity when needed
+        self.r.core.p.power = 0
+        self.r.core.p.powerDensity = 1e9
+        mass = self.r.core.getHMMass()
+        self.r.core.setPowerIfNecessary()
+        self.assertAlmostEqual(self.r.core.p.power, 1e9 * mass)
+
+        # don't take the powerDensity when not needed
+        self.r.core.p.power = 3e9
+        self.r.core.p.powerDensity = 2e9
+        self.r.core.setPowerIfNecessary()
+        self.assertAlmostEqual(self.r.core.p.power, 3e9)
+
 
 class CartesianReactorTests(ReactorTests):
     def setUp(self):