Design of the `MicroXS` and `IndependentOperator` classes

[paulromano]

I'm currently working on refactoring the MicroXS and IndependentOperator classes to support a wider variety of use cases. Namely, the MicroXS class currently stores one-group microscopic cross sections that are to be paired with one-group fluxes. For our FES shutdown dose rate project, we need the ability to do activation calculations using a multigroup flux spectrum / cross sections. Furthermore, the MicroXS and IndependentOperator class are not currently structured to support multiple materials very well (you can have multiple materials in IndependentOperator but only a single MicroXS).

As I'm envisioning a refactor of these two classes, one issue I'm having trouble resolving in my mind is how to handle multiple materials. For example, I can imagine something like the following

# This involves running an OpenMC model to get microscopic cross sections
# by (material, nuclide, reaction, energy group) and fluxes by (material, energy group)
micros, fluxes = get_micros_and_fluxes(model)

op = openmc.deplete.IndependentOperator(materials, micros, fluxes)
timesteps = [...]
source_rates = [...]
integrator = openmc.deplete.PredictorIntegrator(op, timesteps, source_rates=source_rates)

The main issue is: how do we specify the source rates for multiple materials? The options are:

1. The source rate is specified as the actual fixed source rate in [neutrons/sec]. In this case, it would make sense for us to pass fluxes to the independent operator corresponding to a unity source rate, and then all the reaction rates would be scaled accordingly.
2. The source rate is the one-group flux in [neutrons/cm²-s]. This implies that we would have to have a separate instance of the integrator and operator for each material, so the code above would become:

   timesteps = [...]
   for mat, micro_xs, flux in zip(materials, micros, fluxes):
       op = openmc.deplete.IndependentOperator(mat, micro_xs, flux)
       source_rate = flux.sum()
       integrator = openmc.deplete.PredictorIntegrator(op, timesteps, source_rates=source_rate)

Option 2 is arguably closer to how codes like FISPACT work. The multigroup flux provided to the operator has arbitrary units and everything is normalized to the total flux in [n/cm²-s]. Option 1 feels syntactically cleaner and more consistent with CoupledOperator, but it does imply that the units on fluxes matter.

Anyone have thoughts on this? @eepeterson @shimwell @yardasol @pshriwise

[eepeterson]

Thanks for starting this @paulromano! I think if we look at the highest level and go from there the thought process might be something like this:

1. CoupledOperator allows you to deplete an entire model (currently across all materials).
2. Do we want the same interpretation for IndependentOperator? My gut says yes, which means that it too should be able to deplete across all materials (or cells, or mesh voxels, etc) for an entire model.
3. This means the source rate should be for the model i.e. the actual source rate in [neutrons/sec] in IndependentOperator
4. Under the hood, I think we need a helper function to do the Bateman solve for all the cells/materials/voxels independently where each solve uses the local source rate derived from the global source rate (essentially how FISPACT does a single calculation).
5. Maybe in some sense the IndependentOperator becomes a wrapper for other operators. I think it's important to have the functionality to do a FISPACT-like calculation given any arbitrary multigroup spectrum and material composition and then to be able to dispatch that helper function over whatever discretization of your model you want.

I guess when I start thinking about it a bit more I wonder if we need to add an Integrator class to facilitate this or refactor current Integrators to have a method that just does the time integration without needing to know where the burnup matrix came from.

[shimwell]

Being able to easily compare with inventory codes that don't handle transport like would be really helpful. I'm often asked to compared code outputs and was actually thinking of comparing with ALARA as it is open source.

I guess it is double the work but would it be worth fleshing out both approaches as skeleton PRs in case we discover any performance or usability benefits of one approach over the other.

[pshriwise]

Thanks for starting an in-depth conversation about this @paulromano!

My thoughts on this are pretty in-line with @eepeterson's on this. It would be very convenient to have our IndependentOperator class handle the bookkeeping for the domain-specific flux based on a global source term with a lower-level class or function that performs calculations domain-by-domain.

In my head, this would mean that the loop from option 2 above would exist inside the IndependentOperator class after some pre-computation of the source rates for each domain. This would be paired with an internal refactor of IndependentOperator to allow for the depletion/activation of individual domains based on the XS, flux, and a user-specified source rate. Given this, it would be easier for us to enforce the flux units when depleting with the higher-level IndependentOperator functionality since we're handling that ourselves.

As @eepeterson said, an additional step might be to separate burnup matrices from the operators fo for each domain. I'm curious as to whether or not we would be able to "stack" those to perform the Bateman solve more efficiently than repeated calls. In this case, it would be nice to be able to detect when we can avoid duplication of these matrices (i.e. when doing cell-based depletion with cells containing the same material).

A good deal of that is probably me just rephrasing previous thoughts as I understand them 😄 but let me know if I'm missing the mark here.

[paulromano]

Responding to @shimwell:

I guess it is double the work but would it be worth fleshing out both approaches as skeleton PRs in case we discover any performance or usability benefits of one approach over the other.

I don't think it's going to be double the work. It just requires careful thought about the normalization of the multigroup flux being passed to the operator and similarly the source rate passed to the integrator. Clearly in the case where you just have a flux spectrum (absolute units) and you want to activate a material, the source rate will not just be [neutrons/sec] anymore. So, either we need very good documentation and/or some helper functions/classes to ease the burden.

---

Responding to @pshriwise:

In my head, this would mean that the loop from option 2 above would exist inside the IndependentOperator class after some pre-computation of the source rates for each domain.

In fact, IndependentOperator and the integrators are already set up to work with multiple materials, so not too much is needed there. The main change is to allow reaction rates to be calculated with separate flux spectra / cross sections for each material since right now the same ones are used for all materials.

---

Thanks all for the (very prompt) feedback! I agree with the points and it sounds like we're converging around the idea of IndependentOperator handling multiple materials with a global source rate specified for the integrator. I'll start moving in that direction with my refactor; hoping to break it up into smaller chunks as much as possible.

[yardasol]

Hopefully I'm not too late on this @paulromano. To answer your question

how do we specify the source rates for multiple materials?

I think option 1 is more in line with the design goals we initially had for IndependentOperator. I think it's reasonable to assert that units for flux matters! AFAIK that's how the rest of the code behaves as well. Option 2 feels a bit clunky to me, but I'm sure there are use cases that demand it. Maybe as @eepeterson suggested, we could add this functionality through a subclass of IndependentOperator?

In response to the additonal discussion, I have some additional points for consideration in the context of how IndependentOperator was designed:

1. The docstring for IndependentOperator specifies that source-rate normalization option assumes the values being passed are the neutron flux [n/cm^2-s]. This option should probably be changed to make this assumption clearer. The reason I used the term source-rate was that the SourceRateHelper had the correct algorithm set up for calculating the reaction rate by multiplying macroscopic cross sections with a neutron flux. We will probably need to change how we implement the source-rate normalization when adding to multi-group flux and cross section capabilities.

2. It's been a while since I've looked at the code, but I believe IndependentOperator already handles depletion of multiple materials with a single source rate (either a power/power density or a single-group flux). If we look at the User Guide section on Transport-independent depletion,:

   A transport-independent depletion simulation using source-rate normalization will calculate reaction rates for each material independently... A transport-independent depletion simulation using fission-q normalization will sum the fission energy values across all materials into $Q_i$ in Equation (1), and Equation (1) provides the flux we use to calculate the reaction rates in each material.

   Based on this description, it sounds like it wouldn't be super complicated to add functionality where we use a local flux/power/source-rate for different materials, but it's hard to say for sure without analyzing the code itself. Currently, IndependentOperator currently doesn't do anything with the geometry, so adding some machinery to allow localized fluxes/powers/source-rates would be a way to kind of shoehorn the geometry component into the calculation.

3. @pshriwise @eepeterson I think the idea of using a global-source term to create domain-specific fluxes is interesting. Do we need to do a transport simulation to get such quantities? Or is it enough to have the multi-group cross sections and a global source rate?

[paulromano]

I think the idea of using a global-source term to create domain-specific fluxes is interesting. Do we need to do a transport simulation to get such quantities? Or is it enough to have the multi-group cross sections and a global source rate?

A transport simulation would be needed, because you need to know the relative fluxes in each domain based on the source rate. This can easily be streamlined by just extending MicroXS.from_model to support multiple domains (already working on this)

---

I'm planning on breaking down the work as follows:

1. Refactor MicroXS to be based on a numpy array rather than a pandas dataframe. Dataframes are inherently 2D and trying to use them to express data in more than 2 dimensions gets really awkward. Numpy arrays on the other hand allow you to easily express data with N dimensions.
2. Split out MicroXS.from_model into a function that gives you both micro cross sections and fluxes (since we will need the fluxes for IndependentOperator).
3. Modify IndependentOperator to accept fluxes
4. Expand MicroXS and IndependentOperator to support multiple energy groups
5. Expand MicroXS and IndependentOperator to support multiple domains

Question for whoever would like to review my PRs related to this — normally I like to break things up into PRs that are as small as possible, but in this case, successive PRs will keep changing the APIs which might be annoying. Would people prefer that I group all the API-changing behavior into a single PR?

[eepeterson]

Would people prefer that I group all the API-changing behavior into a single PR?

I'd be happy to review this and would say yes to putting it in a single PR. When doing a large refactor it's easier for me to see all the API changes and how they work together in one spot even if it is a big change set.