Impact of RNG on search_for_keff

User Support
1

I’m doing some work regarding search_for_keff and was wondering about the impact of floating point error and random number generation. This post is mostly to get other’s input on expectations for the search/suggestions from their own experiences. My curiosity comes from running the same script, with the same search settings, and seeing slight differences in the search results. For example:
first search

python pincell.py --search --tol 0.8
Iteration: 1; Guess of 2.20e+03 produced a keff of 1.00982 +/- 0.00019
Iteration: 2; Guess of 2.20e+03 produced a keff of 1.00996 +/- 0.00018
Iteration: 3; Guess of 2.19e+03 produced a keff of 1.01141 +/- 0.00019
Iteration: 4; Guess of 2.30e+03 produced a keff of 1.00022 +/- 0.00020
Iteration: 5; Guess of 2.31e+03 produced a keff of 1.00027 +/- 0.00019
Iteration: 6; Guess of 2.29e+03 produced a keff of 1.00167 +/- 0.00018
Iteration: 7; Guess of 2.31e+03 produced a keff of 1.00003 +/- 0.00018
The critical boron concentration achieved with a tolerance of 0.8 was 2310.2362875938898

second search

python pincell.py --search --tol 0.8
Iteration: 1; Guess of 2.20e+03 produced a keff of 1.00982 +/- 0.00019
Iteration: 2; Guess of 2.20e+03 produced a keff of 1.00996 +/- 0.00018
Iteration: 3; Guess of 2.19e+03 produced a keff of 1.01143 +/- 0.00019
Iteration: 4; Guess of 2.30e+03 produced a keff of 1.00075 +/- 0.00018
Iteration: 5; Guess of 2.31e+03 produced a keff of 0.99986 +/- 0.00020
Iteration: 6; Guess of 2.31e+03 produced a keff of 0.99992 +/- 0.00019
Iteration: 7; Guess of 2.31e+03 produced a keff of 1.00055 +/- 0.00020
Iteration: 8; Guess of 2.31e+03 produced a keff of 1.00029 +/- 0.00020
Iteration: 9; Guess of 2.31e+03 produced a keff of 0.99994 +/- 0.00022
The critical boron concentration achieved with a tolerance of 0.8 was 2310.714525935255

third search

python pincell.py --search --tol 0.8
Iteration: 1; Guess of 2.20e+03 produced a keff of 1.00982 +/- 0.00019
Iteration: 2; Guess of 2.20e+03 produced a keff of 1.00996 +/- 0.00018
Iteration: 3; Guess of 2.19e+03 produced a keff of 1.01127 +/- 0.00019
Iteration: 4; Guess of 2.32e+03 produced a keff of 0.99934 +/- 0.00019
Iteration: 5; Guess of 2.31e+03 produced a keff of 0.99991 +/- 0.00020
Iteration: 6; Guess of 2.31e+03 produced a keff of 1.00023 +/- 0.00019
Iteration: 7; Guess of 2.31e+03 produced a keff of 0.99993 +/- 0.00017
The critical boron concentration achieved with a tolerance of 0.8 was 2307.5171242229408

It’s interesting that the first two iterations are identical between searches but then start to have differences on guess 3. In practice, if the final value satisfies the search criterium, then it shouldn’t matter how it gets there, but ideally the fewest number of searches is desired. I think playing with the tolerance would for sure be the way to minimize searches for a given number of particles / batches / etc. Looking at the results, each case is essentially statistically equivalent considering confidence intervals, so I guess these differences don’t matter too much in the end.

I think it might be impossible to generate the same search every time due to changes in the random number stream from running in parallel and/or floating point error causing differences in the guess for each iteration. Just wondering if anyone has thoughts on this (especially which things contribute to differences in searches) and/or recommendations. It seems like erring on the side of more particles will produce searches that end close to each other. I could imagine this would be more susceptible to differences with higher uncertainty in each iteration.

2

I’m also wondering if there’s a way to get the eigenvalue at the critical guess besides re-running the model with that value. It seems the print says

The critical boron concentration achieved with a tolerance of 0.8 was 2307.5171242229408

But it doesn’t report the eigenvalue here. I’ve been running with prints of the solves off, so I’m unsure if it actually does the solve at that concentration or it just stops because of convergence criteria and doesn’t re-run a final critical run (which would be fine).

3

Hi Ligross,

I have played around a little with search function and you are correct that the better your counting statistics (particles/batches/gens) the less variance between the same criticality searches. If you theoretically could run the simulation with counting statistics that result in no variance, then the criticality search would be the exact same between runs as the floating point rounding is almost guaranteed less than your required tolerance or error in k-eff in normal cases. I think in general you want your tolerance to be the same magnitude as your variance you want from your k-eff. You have to be careful here though as tolerance depending on the solve method is relative or absolute, and therefor the units of your criticality parameter can matter when selecting tol (and it might be desirable to manipulate this before running criticality searches). Therefore, the more precision you are looking for in the k-eff the better the tolerance you should use. From what I have gathered from people, if you are looking at keeping a reactor critical for depletion calculations you can get away with less accuracy, but if you’re monitoring things like the reactivity devices within a reactor you will require a lot better statistics. Also, if you’re adjusting something with significant impact on multiple parts of your geometry it can take a long time for the simulations.

As for whether it completes a last transport, I don’t believe it does for the initial guess method as it displays the transports results if not hidden and there is a lot less time for the last guess to be given. It appears to do a type of averaging between the two bounds it found in the final guesses by assuming a line connecting the closest guesses and choosing a point along it closest to the bound that was closest to criticality. If you do the bracket method though I found that it does actually run the last guess, difference between the secant (initial guess) method and the bisect (bracket) method.

Depending on your specific scenario utilizing the python api could be beneficial to speed up the process of searching as you can dynamically input better initial guesses or brackets. For example, if you were completing a boron criticality search following depletion simulations on a reactor with a constantly decreasing reactivity over the lifespan being evaluated, you could improve the bracket method by dynamically adjusting the upper bracket to whatever the previous critical boron concentration was found as you know that reactivity will decrease. If you are really familiar with reactor effects and can expect a maximum reactivity change between each depletion step you could have the lower bracket be in respect to the upper bracket (lower bracket = upper bracket - X) as long as you know that the difference in k-eff for each step will not be greater than X. Lastly, for long simulations where you are completing many criticality searches (such as simulation the entire fuel depletion process for any reactor), it might be useful to run a criticality search function with low counting statistics (1/10th or 1/20th # of particles of your normal particles for example) until it loosely converges to criticality, and then run another criticality search with initial guess being the previously found value and having normal counting statistics which should run for a few extra iterations to converge closer to criticality.

If the reactor’s reactivity does not continuously decrease (CANDU at beginning of cycle for example actually experiences a bump due to plutonium formation) then the initial guess method would be preferred but you should still have the initial guess dynamically change to whatever your previous criticality boron concentration was. There are a few situations that if you think about it you can optimize the searches, which can make a significant difference in long sims that require repeated criticality searches. Lastly, there are options for the bracket method solver, and I would suggest to switch it to brentq unless you want to really evaluate the differences between the solvers.

1 Like
4

According to @pshriwise, using >2 OpenMP threads can lead to differences since floating point addition between threads is not guaranteed to be commutative (A+B =/= B + A). I’ve tested again with 32 MPI processes with one OpenMP thread per MPI process. This is leading to identical searches!

For others wondering how, you can pass the correct run_args to the search like so

run_args = {"mpi_args": ["mpiexec", "-n", "32"], "threads": 1}

which will guarantee the number of processes and OpenMP threads.

1 Like
5

Interesting, I figured it had something to do with the way they use seeds before each run, I would have guessed that the floats had a very high amount of digits being stored (I saw somewhere it’s 8 byte, which should store 16 digits), therefore floating point errors would be insignificant to the variance within the k-eff even in the worst case rounding scenarios.

6

One other thing I’m wondering about is how the two initial guesses are provided – to scipy – for the secant method. It looks like when you don’t provide a bracketed method, openmc says that a secant method is used. It calls scipi.optimize.newton and provides one initial guess. From openmc/search.py

    elif initial_guess is not None:

        # Generate our arguments
        args = {'func': search_function, 'x0': initial_guess}
        if tol is not None:
            args['tol'] = tol

        # Set the root finding method
        root_finder = sopt.newton

    else:
        raise ValueError("Either the 'bracket' or 'initial_guess' parameters "
                         "must be set")

    # Add information to be passed to the searching function
    args['args'] = (target, model_builder, model_args, print_iterations,
                    run_args, guesses, results)

    # Create a new dictionary with the arguments from args and kwargs
    args.update(kwargs)

    # Perform the search
    zero_value = root_finder(**args)

    return zero_value, guesses, results

Looking into the scipy docs, the newton function expects either an fprime (to run Newton) or an x1 that does not equal x0 (to run secant). I can’t figure out from the source code is how it provides the second initial guess – which is a requirement for the secant. If anyone knows how it selects the x1 for secant, please let me know. From my results, it seems like my x0 is the initial guess I give and x1 seems to be close by, but I can’t find the actual code that selects x1

I also wanted to talk more about tolerance and the search method.

I think in general you want your tolerance to be the same magnitude as your variance you want from your k-eff. You have to be careful here though as tolerance depending on the solve method is relative or absolute, and therefor the units of your criticality parameter can matter when selecting tol (and it might be desirable to manipulate this before running criticality searches). Therefore, the more precision you are looking for in the k-eff the better the tolerance you should use.

I’m currently doing searches with tol = 0.5. My current statistics are producing eigenvalues in my search with

Iteration 0 had 2200.0 resulting in 1.0098160667672698+/-0.00018721841437050537
Iteration 1 had 2200.2201 resulting in 1.0099607722920139+/-0.00018356704021688818
Iteration 2 had 2185.0695659388603 resulting in 1.011652646796366+/-0.0001883119546961225
Iteration 3 had 2289.4176258786706 resulting in 1.001903863015256+/-0.00020276848159237997
Iteration 4 had 2309.796005017416 resulting in 1.0000511017611478+/-0.00018601420478047472
Iteration 5 had 2310.358069374611 resulting in 0.9998218612551238+/-0.00020064005290132091
The serach resulted in a criticial guess of 2309.9212990984197

It seems like by your recommendations, my tolerance is a lot higher than the variance (0.5 vs (0.00019)^2), but my search seems okay? I’ve also tried runs with stricter tolerances and they seem to have trouble stopping because the tolerance is too tight. Where does this recommendation derive from exactly? Looking at the docs, scipy.optimize.newton has a tol and rtol parameter, but OpenMC passes tol as part of run args, so I’m thinking my specified tol is absolute.

7

For the way secant selects their next guess, the Scipy website itself states that it just chooses a close guess and checks, but I found this website,

that suggests the method used is basically checking which direction is approaching 0 and using that direction. I have had the secant method fail on me, it has sometimes selected a terrible second guess that was impossible, and I have had iterations go in the wrong direction, not towards criticality, therefore I would warn about the stability of this method over the bracket method, although the bracket method is very slow if you are not really familiar with the problem enough to select tight bounds for initial guesses. I don’t know if there is a way to pass the second guess to the secant method within OpenMC (I haven’t looked for it yet), though I would be interested if there was a way as I think it would add a lot of stability to the process. As for the tolerances, secant I believe uses absolute, but the bracket methods can differ. The suggestion was something I read in one of the forum pages though I don’t remember where. I think it is also a very loose suggestion depending on your scenario. Tolerances are in regards to x, so whatever your parameter you are adjusting is and not the differences in k-eff. So in general, if adjusting your parameter makes a significant change relative to the amount it was adjusted, then you need to be more careful with your tol and vice versa. Example, if 1ppm changed your k-eff by 1mk, versus 1 ppm changing it by 0.01mk, you would need a lot stricter tolerances on the first case as smaller adjustments make greater changes in k-eff. I don’t know an official way to monitor this, I would just always run tests and change the tolerance to fit your specific situation (based on speed and how off your results are from critical), so for your example you show there, you can adjust the parameter by 110 (your entire range), and your k-eff only changes by ~10mk, therefore roughly corresponds with 1 unit/ 0.09mk worth of change over the whole range, or even taking your last two guesses a 1 unit/ 0.41 mk change. Therefore a tol of 0.5 would at most result in a difference of 0.205 mk. So you’re variance in k-eff should be somewhere at this value if you’re wanting to balance computational efficiency. This is all really rough as the effects of a lot of parameters are not linear, and the computational impacts of more iterations or more counting statistics are complicated and coupled but you get the gist of it. For my problem, I deal with boron concentrations of 1.5*10^-05, therefore my tolerances are chosen pretty low (10^-5) (although you could just adjust units here to have a higher tolerance for the secant method, but not for bracket method). Let me know if you have any other questions!

1 Like
8

Ah gotcha, yeah I figured out how mathematically the second guess is generated for the secant method when you don’t supply one and it’s pretty much what you’ve stated. It generates one using a small \epsilon of 1e-4 and then increases (or decreases if negative) the guess by 1+\epsilon

    else:
        # Secant method
        method = "secant"
        if x1 is not None:
            if x1 == x0:
                raise ValueError("x1 and x0 must be different")
            p1 = x1
        else:
            eps = 1e-4
            p1 = x0 * (1 + eps)
            p1 += (eps if p1 >= 0 else -eps)

This might be a fine thing for OpenMC to do, but I still think it would be cool to allow the user to specify their own value. I’ve opened an issue for it and intend to get a PR going next week to add the feature.

1 Like
9

My PR is live, but it looks like another PR involving reactivity control is getting attention (and may supercede my proposal) and could be of interest to you @Jarret