Criticality search method on the Model class

[paulromano]

Description

This PR introduces a new criticality search capability as the Model.keff_search method. This method makes several improvements over the existing search_for_keff function:

• Uncertainty-aware: The new functionality is based on the GRsecant method from a paper by Price and Roskoff (thanks @deanrp2!), which is essentially a secant method but based on an uncertainty-weighted linear regression of recent search points.
• Adaptive batching: Based on the target uncertainty, the algorithm adjusts the number of batches to run. It starts off with few batches and successively ramps up once it gets closer to the target k-eff.
• Simpler user interface: Instead of providing a function that returns a Model class, the keff_search method requires passing a function that makes some change to a mutable object. As the example below shows, it can be as simple as a one-line function.

Fixes #1656

Submitting in draft status to get some community feedback (pinging @church89 @gridley)

Example

Here is a very small example that shows how the new method works:

# Create model of a sphere with pure U235
u = openmc.Material()
u.add_nuclide('U235', 1.0)
u.set_density('g/cm3', 5.0)
sphere = openmc.Sphere(r=1.0, boundary_type='vacuum')
cell = openmc.Cell(fill=u, region=-sphere)
geom = openmc.Geometry([cell])
settings = openmc.Settings(particles=1000, inactive=10, batches=50)
model = openmc.Model(geometry=geom, settings=settings)

# Function that is used during keff search
def modify_radius(x):
    sphere.r = x

result = model.keff_search(modify_radius, x0=0.1, x1=10., output=True, k_tol=1e-3, sigma_final=1e-3)

This produces the following output when I run it:

Iteration 1: batches=40, x=1.00e-01, keff=0.00332+/-0.00004
Iteration 2: batches=40, x=1.00e+01, keff=0.33585+/-0.00123
Iteration 3: batches=20, x=2.98e+01, keff=0.94723+/-0.00615
Iteration 4: batches=20, x=3.02e+01, keff=0.95518+/-0.00524
Iteration 5: batches=20, x=3.06e+01, keff=0.96838+/-0.00526
Iteration 6: batches=20, x=3.16e+01, keff=0.99807+/-0.00412
Iteration 7: batches=303, x=3.17e+01, keff=0.99993+/-0.00124
Iteration 8: batches=519, x=3.17e+01, keff=0.99950+/-0.00095

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 15) 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)

[gridley]

Wow, beautiful. This is the criticality search I've always dreamed of. I have two pieces of feedback on this.

For one, the paper is not open access, so I think it would be better to not refer to Eq. 8 in the docs without saying what it actually is. Maybe some small amount of additional documentation would obviate the need to refer back to the paper?

Secondly, I see no reason why different values for b0 and b1 would be advantageous. It might be good to consider deleting b1 as a parameter.

And lastly, just a question for you: can this be turned on in the course of a depletion calculation? And if so, do you have any thoughts about the steps this would be ran at? Should it be run at every substep? Or perhaps just correctors or predictors for example? There are many ways to implement this.

[paulromano]

@church89 Note that I've just updated the PR to support the use of openmc.lib.

@gridley Thanks for your comments. I've added a little more documentation so it doesn't require a user to go dig through the paper to find that equation. I've also removed b1 as you suggested. As far as use during depletion, that is exactly what @church89 has been working on (via #2693), and in fact this PR is supposed to be in support of that feature. The plan is once this is merged to refactor #2693 to use the new search method here.

[church89]

batches += self.settings.inactive

[church89]

Hi @paulromano, this looks indeed very good and neat. Just a consideration to use this with a depletion calculation and with #2693. In a depletion calculation the openmc.lib model is initialized and therefore the ModelModifier function should be able to access the model stored in memory and modify it there directly. In your example we would need a setattr function to be able to do it or was your idea to let the user crate the function accordingly? In that case modify_radius would have a setattr directly inside.