Refactor MicroXS class and usage in IndependentOperator

[paulromano]

Description

This PR refactors the MicroXS class and corresponding downstream usage in the IndependentOperator class. The full motivation is described in #2580 so I won't go into too much detail here other than to say that now:

• MicroXS no longer subclasses pandas.DataFrame and instead stores a numpy array as the data attribute
• MicroXS can now store multigroup cross sections
• A new get_microxs_and_flux function gives the user a way of getting both (multigroup) microscopic cross sections and corresponding fluxes for a given set of domains.
• The user can now explicitly specify which nuclides/reactions they want instead of always inferring this from the depletion chain
• The IndependentOperator class now requires that you pass it fluxes for each material to be depleted.

Altogether, these changes make it pretty trivial to switch from CoupledOperator to IndependentOperator — just call get_microxs_and_flux to get fluxes and microscopic cross sections and pass them to IndependentOperator with all the other usual info.

Closes #2580

Checklist

• I have performed a self-review of my own code
• I have run clang-format on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[yardasol]

Hi @paulromano,

Overall, I like these changes. See my in-line comments for suggestions.

I initially wrote the below text in an in-line comment for lines 345-348, but decided to put it in my general comments box so we could discuss it without it getting buried in the comment stack.

---

It took me a while to convince myself that this was alright, but it seems that the units work out.

In the initial implementation, _IndependentRateHelper._results_cache was filled by $\sigma^j_i * N_i$ for $j$ reactions and $i$ nuclides. We were able to do this by abusing the structure of the code in OpenMCOperator and ChainFissionHelper. and having the normalization factor output a flux for us in n/cm^2-s based on the equation:

$$ \phi = \frac{P}{\sum_i (Q_i \sigma^f_i N_i)} $$

(https://docs.openmc.org/en/stable/usersguide/depletion.html#equation-fission-q), which would then be multiplied by the cross sections, yielding reaction rates. This allowed us to calculate reaction rates based on just the cross sections, the source rate (power), and the fission Q values.

Now because of the changes from previous PRs removing the need for dilute_initial, this process appears to have some extra steps, with multiplication and division of various factors related to converting barns to cm^2 and vice versa, so I think that the linked equation no longer holds true.

I'm a bit concerned about this, as it means that instead of calculating a flux based on the power, cross sections, fission Q values, and number of atoms, then multiplying that flux by the cross sections to get reaction rates, we are instead relying on MG flux values that we calculate at starting composition to get the reaction rates. These MG flux values stay the same no matter how the composition changes, which isn't going to be accurate for all problems

I've been thinking of how to incorporate the MG cross sections into an equation similar to the one above. I've come up with the following:

$$ \phi_j = \frac{P}{\sum_i (Q_i \sigma^f_{i,j} N_i)} $$

for energy group $j$.

We would calculate $\phi_j$ for all $j$ using _normalization_helper.factor(), then multiply each cross section by it's corresponding $\phi_j$, then sum over all j to get the reaction rates.

While the math does work out, we'd need to have flow control in OpenMCOperator._calculate_reaction_rates(), as IndependentRateHelper.get_material_rates() would return a quantity indexed by reaction, nuclide, AND energy group. We'd also need to modify the ChainFissionHelper to handle the additional dimension in the fission_rates argument.

What are your thoughts on this? Maybe I should open an issue for this?

[yardasol]

See general comments box for big spiel on this

[paulromano]

I'm a bit concerned about this, as it means that instead of calculating a flux based on the power, cross sections, fission Q values, and number of atoms, then multiplying that flux by the cross sections to get reaction rates, we are instead relying on MG flux values that we calculate at starting composition to get the reaction rates. These MG flux values stay the same no matter how the composition changes, which isn't going to be accurate for all problems

If you are using a one-group flux and you have only a single material that is being depleted, the value of the user-provided flux is 100% inconsequential because everything gets normalized in order for the reaction rates to produce the correct amount of power. If you have a multigroup flux, the flux will then determine the proper weighting between the provided multigroup cross sections. Similarly, if you have multiple materials being depleted in IndependentOperator, the per-material fluxes provide proper weighting between the different materials. Without providing fluxes and using them when calculating the reaction rates, there's no way we can get the proper weighting. Note that in all cases the power is still used to adjust the reaction rates at each step.

You are correct that perhaps our equations need updating -- I'll look at what we have written and think about that further.

While the math does work out, we'd need to have flow control in OpenMCOperator._calculate_reaction_rates(), as IndependentRateHelper.get_material_rates() would return a quantity indexed by reaction, nuclide, AND energy group. We'd also need to modify the ChainFissionHelper to handle the additional dimension in the fission_rates argument.

I'm not sure I understand where you're going with this. Maybe we should hash it out over a call 😄

[yardasol]

@paulromano

If you are using a one-group flux and you have only a single material that is being depleted, the value of the user-provided flux is 100% inconsequential because everything gets normalized in order for the reaction rates to produce the correct amount of power. If you have a multigroup flux, the flux will then determine the proper weighting between the provided multigroup cross sections.

This makes sense, but wouldn't the weighting change over time as the neutron spectrum shifts?

[paulromano]

@yardasol Yes, if the material composition changes enough that it perturbs the neutron flux, the weight can change. That's one of the reasons we need to run transport solves at every step for fission reactor depletion. However, for fusion systems, in a lot of cases the material is getting irradiated and activated but doesn't change enough to actually perturb the flux (effectively, there's no "burnup"). In those cases the original weighting from time t=0 should hold just fine.

[eepeterson]

As always, really great stuff @paulromano! I have a few questions/comments as I'm digging into this before a full review:

1. Since we're doing some refactoring now, does it make sense to "promote" MicroXS to outside the depletion module? I see a few reasons for this:

• The concept/use of microscopic cross sections is more general than just for depletion.
• It may help for future refactoring of the mgxs module.
• It would be nice for helping with plotting/understanding cross sections generally, especially what goes into generating multigroup cross sections and the differences between fission and fusion data needs.

If the answer to 1. is "yes", then I think it also makes sense to promote the depletion dependencies within microxs.py: namely the chain file handling in chain.py (and it's dependency on nuclide.py) as well as the _find_cross_sections and _get_nuclides_with_data methods within coupled_operator.py by moving them to the Model class and data/library.py? respectively. I think this helps minimize code in the depletion module prior to migrating it to C++ and provides a more organized structure.

2. Can we add a couple "TODO"s to the MicroXS class? One to be able to generate a MicroXS object directly from cross section libraries rather than running a transport simulation and another to be able to generate a burnup matrix directly from the MicroXS object.

3. I foresee some memory problems for large R2S problems in the way get_microxs_and_flux works, so it might be worth thinking about how to tally the non-zero burnup matrix entries directly for every domain as well as the total flux in each domain.

Interested to hear feedback on these points and how much of it makes sense to add in this PR vs another PR vs never.

[paulromano]

1. does it make sense to "promote" MicroXS to outside the depletion module?

Yes, I think there's probably a good argument for this. I'd prefer to keep any sort of file restructuring for another PR just to limit the scope of each of them. Is that OK with you?

2. Can we add a couple "TODO"s to the MicroXS class? One to be able to generate a MicroXS object directly from cross section libraries rather than running a transport simulation and another to be able to generate a burnup matrix directly from the MicroXS object.

Already have plans for this 😄 But yes, I'll add some TODOs in this PR for good measure

3. I foresee some memory problems for large R2S problems in the way get_microxs_and_flux works, so it might be worth thinking about how to tally the non-zero burnup matrix entries directly for every domain as well as the total flux in each domain.

You'll have to elaborate more, but I guess the only opportunity I see is if you were to do each nuclide separately since we know that some reactions won't apply to all nuclides. Is that what you were hinting at?

[eepeterson]

Sounds great to me - fine to leave the restructuring for another PR and glad to hear you're already way ahead of me on more features for MicroXS.

As for my point 3 re: memory issues: Yes, essentially doing each nuclide separately, but through some combination of preprocessing continuous energy cross sections into continuous energy transmutation fractions or doing it on the fly through something like a "reaction sum" score since the transmutation pathways should be fewer than the total possible reaction pathways.

[eepeterson]

Thanks for tackling this @paulromano I'm happy with merging this if @yardasol is.

[paulromano]

@yardasol @shimwell any further comments on this PR?

[shimwell]

No

@yardasol @shimwell any further comments on this PR?

Looks great to me LGTM

[shimwell]

I think the small conflict on this PR might be due to me replacing OrderedDicts with {}, sorry about that

[yardasol]

@paulromano apologies for the late response! I was in the middle of moving apartments when you pinged me. Have you come up with an equation to replace equation (1) in our doc pages? It would be good if we could have this before we merge.

[paulromano]

@yardasol Good suggestion -- I've just updated that equation. Let me know what you think.

[yardasol]

@paulromano I like your update but we should perhaps state how the normalization factor is used to get reaction rates.

[paulromano]

@yardasol I just made an update to clarify how the normalization factor is used. Hopefully that's clearer now.

[yardasol]

@paulromano This is coming together. I think you are missing some terms in the new equation (specifically the number / (1e24 * volume) term that is called atom_per_bcm.

Also, I feel this may have been covered at this point in one of the previous PRs, but why are the units of flux in the MicroXS class n-cm/src? I recall there being a different term for this quantity but it is proportional to the flux by a constant factor.

[eepeterson]

Also, I feel this may have been covered at this point in one of the previous PRs, but why are the units of flux in the MicroXS class n-cm/src? I recall there being a different term for this quantity but it is proportional to the flux by a constant factor.

This unit of flux is just different from the "normal" flux units in that it is not yet divided by the volume over which it was integrated (tallied). In this particular use case it is because the OpenMCOperator requires the get_material_rates method to return reaction rates in units of [(reactions/src)*b-cm/atom] since it itself divides by the volume within the _calculate_reaction_rates method.

[paulromano]

I think you are missing some terms in the new equation (specifically the number / (1e24 * volume) term that is called atom_per_bcm.

@yardasol The number density is already accounted for in the $N_{i,m}$ term (total number of atoms). When we tally "reaction rates" now, because we are using multiply_density=False, the tally values do not account for the number density of each nuclide (since it is $\phi\sigma$, not $\phi\Sigma$). Multiplying by atom_per_bcm in OpenMCOperator._calculate_reaction_rates gives us the reaction rate accounting for the density of each nuclide. Does that make sense?

[yardasol]

@eepeterson @paulromano thanks for those explanations. I was just asking because these units might make it more difficult for someone who wanted to manually provide those fluxes (i.e. from a separate code, which the original implementation of MicroXS supported),

[yardasol]

Good work @paulromano, just one final comment, but otherwise I think this is good to go. Looks like we are having some RTD issues, so should I hold off on merging?

[yardasol]

Maybe add a note specifying that the units of the flux are not the expected [n/cm^2-s]

[paulromano]

I'm hesitant to make a note about the units here -- if we were to change it in the class, it will be easy to forget to change it here too. The IndependentOperator class docstring is explicit on what units are expected.

[paulromano]

Looks like we are having some RTD issues, so should I hold off on merging?

The RTD issue was fixed in #2667; I merged the develop branch into this one and it looks like RTD is happy now. @eepeterson @yardasol if you guys are happy, feel free to hit the big green button!

[eepeterson]

all checks passing and a couple approvals - seems good to me!

[Ebiwonjumi]

Hi @paulromano, concerning the get_microxs_flux function, it will be good to save the reaction rates and flux values from the temporary run to a statepoint file for restart purposes.

[eepeterson]

@Ebiwonjumi thanks for chiming in! This is an already merged PR so if there is additional functionality you'd like to see could you create a new issue if it isn't already there?

[Ebiwonjumi]

Got it. Thanks @eepeterson