Add reactivity control to coupled transport-depletion analyses

[church89]

Description

This pr adds the capability to perform batch-wise reactivity control during a coupled transport-depletion analysis.
It effectively runs an iterative search_for_keff sub-step algorithm to maintain keff equal to 1 (or some other user defined targets) as a function of some user defined model parameter to update, before calculating the BOS from the transport run and proceed with the depletion step.
The model variable to parametrize during the sub-step can be of geometrical (e.g. control rod position) or material (e.g. refuel) nature, acting on a openmc.Cell and openmc.Material, respectively.

A batch-wise scheme can be added to an integrator instance and a coupled transport-depletion simulation can be run like this:

...
...
op = openmc.deplete.CoupledOperator(model)
integrator = openmc.deplete.CECMIntegrator(op, timesteps , power)
integrator.add_batchwise(cell, 'translation', axis = 2, bracket = [-4, 4], bracket_limit = [-100,20])
integrator.integrate()

where cell is an openmc.Cell conveniently defined that can be translated along the z-axis during the run.

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)

[drewejohnson]

@church89 thank you for this lengthy and useful PR! I think this is going to unlock a lot of doors for people and provide a lot of value. I haven't done a full review but I wanted to get some of my larger thoughts out first, make sure we understand the vision and goal, and then refine my review as we go.

Some simple things

• Please re-evaluate the docstrings you've added to align with the numpydoc spec Otherwise our built documentation isn't going to process this properly
• Please run a python formatter like black to make sure your code adheres to pep8 formatting style like the dev guide requests
• Please consolidate property setters / getters / validators. I think it's going to be easier for us in the long run to maintain / add validators if we do validation on things like self.tol if all validation is done in property setters

Maybe a larger lift: I think the high-level class and method name, using batch wise, might not communicate what we want to communicate. My rationale is I believe this engages after a full transport simulation, while openmc uses the "batch" terminology to mean generations of particles in a given transport simulation - docs

so users who already are configuring an openmc run to have so many inactive / active generations per batch may expect this search to be performed more often than it is actually performed.

But, based on your PR message

It effectively runs an iterative search_for_keff sub-step algorithm to maintain keff equal to 1 (or some other user defined targets) as a function of some user defined model parameter to update, before calculating the BOS from the transport run

and my non-exhaustive review, I could be wrong. If we are making changes to the model during batches of active generations, then this name makes sense. But if we are not, I would encourage us to consider new names.

I hope to have time in the next week or so for a more detailed review, but I wanted to get something out the door to let you know I've seen this, I believe it's valuable, and I have some (hopefully) low-hanging requests

[drewejohnson]

Oof a 2% relative tolerance is high. But we've got a stochastic code, and probably wanting to get the simulation to converge reasonably quickly.

I'd like to see a tighter tolerance but I understand that might not be tenable.

[church89]

Hi @drewejohnson, thanks for starting this review. I've begun to address your comments/suggestions.

I can see how batchwise might create some confusion. The routine acts in between a transport and a depletion step, so an alternatives could be stepwise to maintain the general purpose.
Otherwise, we might think of linking the concepts of depletion and criticality, as ultimately is what it does, like depletion_critical , criticality_depletion or something similar.

[drewejohnson]

Thank you for incorporating some of the changes. I want to make one last call for the following items to be either incorporated or rejected with explanation.

Infinite loop

I've personally been bitten by an unterminated while loop and it can be deeply frustrating. And, since users of this feature are likely going to be running expensive simulations, there is a likelihood that someone expends some unnecessary computational resources because there is no hard upper limit on iterations.

If there is some other way that I'm not seeing where the main loop will terminate if it cannot find a satisfactory solution, please let me know.

reactivity controller setter

The add_reactivity_control method is nice, I'm not trying to deny that. But it does not make it easy for people to develop their own reactivity controller outside of openmc and then use that in a simulation. They would need to modify the source code (either locally or through a PR here) to add their specific flavor of controller. Or they push an issue here and wait for someone else to make those changes.

If the reactivity_controller is a publicly settable attribute, then people have the freedom to use their own controller without needing to make changes to the core of openmc. It can (and should) be used by the add_reactivity_control helper method.

search.py

These changes don't feel intrinsically related to the PR at hand. They feel valuable, but the relationship is not clear to me. There is some issue you are trying to workaround that isn't transparent to me.

I also have some concerns about variable return signature (two or three arguments) and would welcome other opinions on the matter.

[church89]

@drewejohnson
Thank you for your comments. Regarding your concerns:

• Infinite loop: As long as the search_for_keff bracketed algorithm returns two arguments (meaning that the search have failed to bound the root) the adaptive algorithm will try to move the bracket range to bound the target. When that happens the search_for_keff returns three arguments and the first exit point to the while loop.
  If the adaptive algorithm moves the bracket range outside of the user defined limits (e.g. the entire span of a control rod), the closest limit is set as the root and the while loop exits as well.
  The only case the loop won't exit easily is if the adaptive algorithm won't bound the target but still remain inside the user defined limits, this might take awhile depending on user settings and simulation model but eventually should make it.
  Inside the loop there is also another check condition that is: when search_for_keff returns two arguments, (guesses and keffs) and keffs difference is below 1 pcm, the algorithm is already very close to target and the root will be set as the closest guess to target, exiting the loop. This 1 pcm can be user defined, but one should not go lower than 1, otherwise there is a risk of overshooting the adaptive bracket out of user limits and invalid the method

• Reactivity control setter: Generally I like your idea and I think it makes sense, in practice I think it is not so straightforward and by letting users define their own controllers there is still a high risk of not being able to prevent misleading code. We would need to think of some very general helper methods that would need to include any user ideas and that might not be an easy task. What cases were you thinking that are not already covered by existing controllers?

• Search method: These changes here are related to this PR. As explained above in the infinite loop, the search_for_keff must be able to fail (returning two arguments instead of three) to enable the adaptive algorithm and find the right bracket range the bounds the root. The run_in_memory argument instead is simply to enable search_for_keff to work with a model instantiated in memory. This is also related, as the operator class runs the model in memory.

[church89]

@drewejohnson
I thought a bit more about the reactivity control definition and added a setter method that hopefully goes in the right direction.
Now users can instantiate directly their own ReactivityController class similar to what happens for the TransportOperator ones.
So now a user can either set directly aReactivityController based on the current implementation of the reactivity_control.pymodule or defined its own.
In the setter method there are a few checks for inheritance and attributes, but there is no check for the different methods needed to build a reactivity controller (and here is where to tricky part might come).
Just would like to have your feedback before moving on and if that is pretty much what you were thinking.
Thanks!

[drewejohnson]

Apologies again for my delay. Life has been a bit hectic and I finally found time to review this.

Thank you for adding the reactivity control setter. I think it enforces too many checks and doesn't allow users to expand beyond what OpenMC can do already without source code changes. Please consider requiring the setter to only check that the object is an instance of the base ReactivityController. If there are concerns about allowing the users to construct those subclasses properly, what checks and validation steps can be added to the subclasses? e.g., to ensure that CellReactivityController.cell is always a Cell subclass?

Regarding the infinite loop,

The only case the loop won't exit easily is if the adaptive algorithm won't bound the target but still remain inside the user defined limits

This is not an uncommon situation. Consider control rod movement through depletion. It is conceivable that, towards the end of life, there is no way to remove the control rods enough from the reactor to achieve a critical state with sufficiently spent fuel. Same for other reactivity control mechanisms where a poison is modified (control rods, control drums, soluble boron, etc.)

I am going to signal for an additional reviewer to help close the gap on the infinite loop and API changes to openmc/search.py.

[drewejohnson]

Why only these four control mechanisms? What about the interplay between the Operator and the ReactivityController requires only these four mechanisms?

For light water reactors, adding and changing soluble boron concentration might be akin to "refueling" but we're not adding fuel there so someone could set attrib_name = "dillution" (for example).

[church89]

Hi @drewejohnson and thanks again.

Apologies again for my delay. Life has been a bit hectic and I finally found time to review this.

Thank you for adding the reactivity control setter. I think it enforces too many checks and doesn't allow users to expand beyond what OpenMC can do already without source code changes. Please consider requiring the setter to only check that the object is an instance of the base ReactivityController. If there are concerns about allowing the users to construct those subclasses properly, what checks and validation steps can be added to the subclasses? e.g., to ensure that CellReactivityController.cell is always a Cell subclass?

I've reduced the reactivity controller classes to only two: CellReactivityController and MaterialReactivityController, which both inherit from the abstract ReactivityController, and update the setter. It should be possible in principle now to have the users defined their own controllers.

Regarding the infinite loop,

The only case the loop won't exit easily is if the adaptive algorithm won't bound the target but still remain inside the user defined limits

This is not an uncommon situation. Consider control rod movement through depletion. It is conceivable that, towards the end of life, there is no way to remove the control rods enough from the reactor to achieve a critical state with sufficiently spent fuel. Same for other reactivity control mechanisms where a poison is modified (control rods, control drums, soluble boron, etc.)

That's exactly the reason we have to set a bracket_limit: if one of the bounds (upper or lower) gets hit the while loop exits and the search is stopped.
In case of a removable control rod, you would set this bracket_limit to be the active length of your rod inside the core, so that when the it gets pulled out completely the routine knows that there is not margin anymore and needs to stop.

[paulromano]

Finally getting around to giving this a proper review. Thanks for the work on this PR and for your patience @church89. And also many thanks to @drewejohnson for providing very thoughtful comments that have helped shaped the PR thus far.

Overall I have a few thoughts:

1. I personally like the idea of giving users more control over the criticality search and keeping it general rather than having a bunch of subclasses of a ReactivityController class. Right now, if a user wants to do some other search, they are forced to write a class that implements the ReactivityController interface. I feel like this could be much simpler overall. There could be a generic reactivity controller that accepts a function that takes one argument (the search value) and makes a corresonding change in the model via openmc.lib. So for example, if the user wanted to translate cell 50 along the y axis, the function could be as simple as:

   def translate_cell(value):
       cell = openmc.lib.cells[50]
       cell.translation = (0, value, 0)

   This may arguably be less fool-proof than providing specific controllers but we can try to make things more obvious for users by providing examples (in the docstrings and/or example notebooks).
2. I agree with @drewejohnson's earlier comment about the changes in search_for_keff. Let's separate those out into a new PR to get them merged sooner.
3. I've been wanting to make some other changes in search_for_keff, namely to use the scipy.optimize.root_scalar function that can dispatch to brentq, bisect, etc. This would make it future-proof (if a new method is added, we don't need to update OpenMC). It probably makes sense to get that in before we merge this PR since it might affect the interface, so I'll try to get that PR submitted soon.

[gridley]

On the use of root_scalar: have you all experimented with using that in a function with noise? I just wanted to share this paper since it presents a method for adjusting some parameter (drum angle in this case) to find the root of a function you can only get samples of with noise. It's not equivalent to applying secant method or bisection to the samples.

So, that being the case, I would argue that using the root_scalar API would not be the move here necessarily.

https://www.researchgate.net/publication/372814209_Method_for_control_drum_position_critical_search_with_Monte_Carlo_codes

And lastly, I agree that with so much robust C API functionality already in place, parametrizing the reactivity search around a scalar passed to C API calls implemented in a python callback would be the cleanest and most general way to implement this!

[church89]

Thanks @paulromano for getting back to this PR. Looks reasonable to me to separate the search_for_keff bit and then maybe rename the ReactivityController to something like OperatorReactivityController or DepletionReactivityController that wraps around it.
Let me know if there is anything I can do to help you with the search_for_keff PR.