reviewer suggestion on m and lfs

Author: shimwell

openmc/deplete/chain.py

@@ -134,10 +134,11 @@ def _read_activation_products(evaluation):
     Returns
     -------
     dict
-        ``{mt: {(Z, A, lfs): yield_func}}`` where yield_func is a
+        ``{mt: {(Z, A, m): yield_func}}`` where yield_func is a
         :class:`openmc.data.Tabulated1D` giving the yield as a function of
-        energy. lfs is the level number flag (0=ground, 1=first metastable,
-        etc.)
+        energy. m is the isomeric state number (0=ground, 1=first metastable,
+        etc.), derived by sequentially numbering the excited states in order
+        of increasing LFS.
     """
     results = {}
 
@@ -207,7 +208,22 @@ def _read_activation_products(evaluation):
                         openmc.data.Tabulated1D(energy, yield_)
 
         if products:
-            results[mt] = products
+            # Remap LFS level numbers to sequential isomeric state numbers.
+            # LFS in MF9/MF10 refers to the energy level, not the isomeric
+            # state number.  Isomeric states are numbered sequentially:
+            # ground (LFS=0) -> m=0, then m=1, m=2, ... for excited states
+            # in order of increasing LFS.
+            za_groups: dict[tuple[int, int], list[int]] = {}
+            for (Z, A, lfs) in products:
+                za_groups.setdefault((Z, A), []).append(lfs)
+
+            remapped = {}
+            for (Z, A), lfs_values in za_groups.items():
+                lfs_values.sort()
+                for m, lfs in enumerate(lfs_values):
+                    remapped[(Z, A, m)] = products[(Z, A, lfs)]
+
+            results[mt] = remapped
 
     return results
 
@@ -230,7 +246,7 @@ def _add_isomeric_reactions(nuclide, rx_name, q_value, iso_products,
     q_value : float
         Q value of the reaction in [eV]
     iso_products : dict
-        ``{(Z, A, lfs): yield_func}`` from :func:`_read_activation_products`
+        ``{(Z, A, m): yield_func}`` from :func:`_read_activation_products`
     decay_data : dict
         Decay data keyed by nuclide name
     missing_rx_product : list
@@ -241,8 +257,8 @@ def _add_isomeric_reactions(nuclide, rx_name, q_value, iso_products,
     # Map product states to daughter names, filtering missing nuclides
     yield_funcs = {}  # {daughter_name: yield_func}
     ground_state = None
-    for (Z, A, lfs), yield_func in iso_products.items():
-        daughter = gnds_name(Z, A, lfs)
+    for (Z, A, m), yield_func in iso_products.items():
+        daughter = gnds_name(Z, A, m)
 
         if daughter not in decay_data:
             orig_daughter = daughter
@@ -251,15 +267,15 @@ def _add_isomeric_reactions(nuclide, rx_name, q_value, iso_products,
                 missing_rx_product.append((parent, rx_name, orig_daughter))
                 continue
 
-        # If multiple LFS map to same daughter, keep the one with
+        # If multiple states map to same daughter, keep the one with
         # larger max yield (shouldn't normally happen)
         if daughter in yield_funcs:
             if np.max(yield_func.y) > np.max(yield_funcs[daughter].y):
                 yield_funcs[daughter] = yield_func
         else:
             yield_funcs[daughter] = yield_func
 
-        if lfs == 0:
+        if m == 0:
             ground_state = daughter
 
     if not yield_funcs:

tests/unit_tests/test_deplete_chain.py

@@ -598,14 +598,14 @@ def test_reduce(gnd_simple_chain, endf_chain):
     <isomeric_branching reaction="(n,gamma)">
       <energies>1e-5 0.0253 1e3 1e5 1e6 1e7</energies>
       <product nuclide="Am242">0.9 0.9 0.87 0.84 0.74 0.52</product>
-      <product nuclide="Am242_m2">0.1 0.1 0.13 0.16 0.26 0.48</product>
+      <product nuclide="Am242_m1">0.1 0.1 0.13 0.16 0.26 0.48</product>
     </isomeric_branching>
   </nuclide>
   <nuclide name="Am242" half_life="57672.0" decay_modes="1"
            decay_energy="0.0" reactions="0">
     <decay type="beta-" target="Cm242" branching_ratio="1.0"/>
   </nuclide>
-  <nuclide name="Am242_m2" half_life="4449312000.0" decay_modes="1"
+  <nuclide name="Am242_m1" half_life="4449312000.0" decay_modes="1"
            decay_energy="0.0" reactions="0">
     <decay type="IT" target="Am242" branching_ratio="1.0"/>
   </nuclide>
@@ -630,13 +630,13 @@ def test_isomeric_branching_chain_roundtrip():
         energies, products = am241.isomeric_branching['(n,gamma)']
         assert len(energies) == 6
         assert 'Am242' in products
-        assert 'Am242_m2' in products
+        assert 'Am242_m1' in products
 
         # Check energy-dependent values
         assert np.isclose(products['Am242'][0], 0.9)
         assert np.isclose(products['Am242'][-1], 0.52)
-        assert np.isclose(products['Am242_m2'][0], 0.1)
-        assert np.isclose(products['Am242_m2'][-1], 0.48)
+        assert np.isclose(products['Am242_m1'][0], 0.1)
+        assert np.isclose(products['Am242_m1'][-1], 0.48)
 
         # Write and re-read
         chain.export_to_xml('chain_iso_rt.xml')

tests/unit_tests/test_deplete_nuclide.py

@@ -360,7 +360,7 @@ def test_isomeric_branching_xml_roundtrip():
     energies = np.array([1e-5, 0.0253, 1e3, 1e5, 1e6, 1e7])
     products = {
         "Am242": np.array([0.9, 0.9, 0.87, 0.84, 0.74, 0.52]),
-        "Am242_m2": np.array([0.1, 0.1, 0.13, 0.16, 0.26, 0.48]),
+        "Am242_m1": np.array([0.1, 0.1, 0.13, 0.16, 0.26, 0.48]),
     }
     nuc.isomeric_branching["(n,gamma)"] = (energies, products)
 
@@ -374,7 +374,7 @@ def test_isomeric_branching_xml_roundtrip():
     assert iso_elems[0].find("energies") is not None
     product_elems = iso_elems[0].findall("product")
     assert len(product_elems) == 2
-    assert {p.get("nuclide") for p in product_elems} == {"Am242", "Am242_m2"}
+    assert {p.get("nuclide") for p in product_elems} == {"Am242", "Am242_m1"}
 
     # Read back from XML
     nuc2 = nuclide.Nuclide.from_xml(elem)