Introduce the alpha-k iteration

[ilhamv]

Description

This PR adds an alpha-eigenvalue (aka time-eigenvalue) transport capability. It modifies the existing k-eigenvalue power iteration to converge to the fundamental (right-most, algebraically largest) alpha mode.

• The k-eigenvalue quantity is still tallied, but it now converges to unity.
• The capability can be turned on by setting settings.alpha_mode = True.
• It works on both CE and MG modes.
• The PR builds up on top of the implementation in [Variansyah 2020].
• settings.prompt_only mode toggle is added, which is useful for getting the prompt alpha mode in subcritical systems.
• Besides the alpha-effective, other related quantities of interest are also reported: multiplication factor (the k-effective weighted by the fundamental alpha-eigenfunction), delayed fraction, and mean lifetime.
• The feature is verified against:
  • Variations of the Godiva sphere [Cullen 2003],
  • Subcritical and delayed-supercritical two-group infinite homogeneous media [Zoia 2015], and
  • 361-group infinite homogeneous media representing infinite PWR pin-cells.
• The verification test problems are available here.

Checklist

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

[paulromano]

Hi @ilhamv. Thanks for submitting this pull request. I apologize I haven't had time to take a thorough look through it yet, but I'm hoping to get to it by the end of the week.

[paulromano]

@ilhamv I really apologize it's taking me so long to get you some feedback on this draft PR. Unfortunately it's been hard to find the block of time to sit down and try to fully digest this, especially since I'm only somewhat familiar with alpha-eigenvalue method and need to read through the literature as well. In any event, I have started going through the code and overall it looks pretty good from what I can tell. A few pieces of general feedback:

• It would be very helpful if in the implementation you could put references to equation numbers from published papers in the comments. For example, it took me a while to figure out that the weight modification that happens in Particle::event_advance is an alternative to having a time-absorption cross section. Once I stumbled upon this paper (not sure if that is the original reference for the weight modification), it became immediately obvious what was going on. Before I try to decode every part of the algorithm, you might be able to save me some time since you know where everything comes from :)
• Some of the data that is stored per nuclide, per precursor group (e.g., alpha_lambda) could be stored in a 2D array rather than a vector of vectors. We use the xtensor package elsewhere in the code for multidimensional arrays. For a 2D array, you'd have a xt::xtensor<double, 2>.
• You'll need to add one or more tests that utilize the alpha eigenvalue calculation so that we can ensure that it doesn't become broken over time. Have a look at some of the tests in tests/regression_tests and let me know if you have specific questions.
• Looks like there is now a conflict with develop because of some recently merged PRs. You'll need to either merge the develop branch into yours or rebase this branch (whichever you feel more comfortable with).

[nelsonag]

Morning @ilhamv, this looks like an interesting PR! I just took a first glance and had the following. I can do a deeper dive after the code is merged, etc.

• The documentation rst files should be updated in this or an immediate follow-on PR.
• I am not a fan of the settings.mg_speeds parameter. This is because either (1) that information could be sufficiently derived from the group structure (if an estimate of the average speed within a group is sufficient), or (2) we already have a temperature- and group-dependent inverse-velocity parameter in the MGXS library which you should use instead.

[LiZhulun]

@ilhamv I really apologize it's taking me so long to get you some feedback on this draft PR. Unfortunately it's been hard to find the block of time to sit down and try to fully digest this, especially since I'm only somewhat familiar with alpha-eigenvalue method and need to read through the literature as well. In any event, I have started going through the code and overall it looks pretty good from what I can tell. A few pieces of general feedback:

• It would be very helpful if in the implementation you could put references to equation numbers from published papers in the comments. For example, it took me a while to figure out that the weight modification that happens in Particle::event_advance is an alternative to having a time-absorption cross section. Once I stumbled upon this paper (not sure if that is the original reference for the weight modification), it became immediately obvious what was going on. Before I try to decode every part of the algorithm, you might be able to save me some time since you know where everything comes from :)
• Some of the data that is stored per nuclide, per precursor group (e.g., alpha_lambda) could be stored in a 2D array rather than a vector of vectors. We use the xtensor package elsewhere in the code for multidimensional arrays. For a 2D array, you'd have a xt::xtensor<double, 2>.
• You'll need to add one or more tests that utilize the alpha eigenvalue calculation so that we can ensure that it doesn't become broken over time. Have a look at some of the tests in tests/regression_tests and let me know if you have specific questions.
• Looks like there is now a conflict with develop because of some recently merged PRs. You'll need to either merge the develop branch into yours or rebase this branch (whichever you feel more comfortable with).

Hi @paulromano .I did some sub-critical benchmark tests with this code, and it took me a long time. During this time I discovered a very important problem. Only the Godiva and Jezebel benchmarks mentioned by Ilham can run stably. There is a "segment fault(core dump)" when running core and component benchmarks (even OpenMC's own Lattice). When debugging with the GDB tool, bug.txt appears that the program is having a problem when it calls "Particle::event_advance()". But as a C++ beginner, I don't see any problems. Can you give me some suggestions for modification (or debugging)?

[ilhamv]

Hi @LiZhulun. Thank you for pointing a bug! If you don't mind, can you please send me (email: [email protected]) your input files?

Also @paulromano and @nelsonag, I'm so sorry for my super delayed response on this PR... Thank you very much for the feedback; I should be able to work on this sometime this month.

[nelsonag]

Hi @ilhamv, I'm just curious the progress of this PR. Thanks!

[ilhamv]

Hi @nelsonag, I'm so sorry for the very slow progress; I still want to work on this PR! I'm targetting to get some progress going starting next month. Is there any particular concern/suggestion that I need to be aware of with regard to the more recent OpenMC dev activities? Thanks!

[ilhamv]

Hi, @paulromano and @nelsonag

I rebased, resolved conflicts, and force-pushed the branch. (I hope it is done correctly...)

I applied @nelsonag suggestion on using MGXS-inverse_speed instead of mg_speeds. It's interesting to realize that multigroup-neutrons may move with different speeds in different materials!

Several items that I'm still working on for the feature:

1. Add more descriptive comments, with references to related papers.
2. Use xtensor instead of vector of vector. I used vector of vector because I allowed different numbers of precursor groups for fissionable nuclides (different sizes for the 2nd dimension). Now I'll assume that all fissionable nuclides have the same number of precursor groups.
3. Regression tests.
4. Add collision and absorption estimators for global alpha-eigenvalue calculation (similar to those in k-eigenvalue).
5. Standardize alpha-eigenvalue solution printouts and hdf5 datasets (to those of k-).
6. Replace continuous weight capture (CWC) with implicit capture. Continuously changing particle weight complicates track-length mesh-tally mechanisms; and compared to usual implicit capture, CWC is more efficient ONLY if the medium is highly absorbing. (7) Documentation rst files?

[ilhamv]

I was trying to write a regression_test for the feature and found an issue: Originally, OpenMC uses the running average of k as the keff for the simulation in the next generation. I see that this helps in improving convergence. However, I found that this may interfere in alpha-k iteration, such that the k and alpha may not converge to the desired values (1 and the actual alpha, respectively). So, I changed the code so that we use the latest k_generation--instead of the running average--as the next keff (which is the normal way?), and the alpha-k iteration works as expected. Unfortunately, as one would expect, this changes the results and fails many of the regression tests.

What do you think, @paulromano @nelsonag? Thank you!

[ilhamv]

Hi @paulromano and @nelsonag,

I think the PR is ready for review.

Summary

It has been a while since the PR was first opened. So, here is the summary:

• The PR adds an alpha-eigenvalue (aka time-eigenvalue) transport feature. It modifies the existing k-eigenvalue power iteration to converge to the fundamental alpha mode.
• The "k-eigenvalue" quantity is still tallied, but it now converges to unity.
• Users just need to add settings.alpha_mode = True on their k-eigenvalue input decks to turn on the feature.
• The feature works on both CE and MG modes.
• The PR builds up on top of the implementation in [Variansyah 2020].
• settings.prompt_only mode toggle is added, which is useful for getting the prompt alpha mode in subcritical systems.
• Besides the alpha-eff, other related quantities of interest are also reported: multiplication factor (the "k-eff" weighted by the fundamental alpha-eigenfunction), delayed fraction, and mean lifetime.

Test suite

The feature is tested on:

• Variations of the Godiva sphere [Cullen 2003],
• Subcritical and delayed-supercritical two-group infinite homogeneous media [Zoia 2015], and
• 361-group infinite homogeneous media representing infinite PWR pin-cells. The test problems are available here.

Key modifications

Key modifications include:

• Adding Particle.speed() for MG mode,
• Continuous weight correction during particle event_advance,
• Redefining nu_fission in tallying keff (track length, collision, and absorption),
• Time-corrected delayed fission neutron emission at fission site creation,
• Alpha eigenvalue update procedure (in eigenvalue.cpp), which requires integral tallies of density (inverse speed), prompt production, and nuclide- and precursor-wise delayed productions.
• Alpha-source site creation and storing into fission bank if fission production < time source (-alpha/v). In this case (simulation::store_alpha_source = true), fission sites are stored in the fission bank instead. In other words, we effectively switch from the iterative scheme in [Variansyah 2020] to [Shim 2015]. Note that the switch is done on the fly and typically occurs in simulations of prompt alpha mode of deeply subcritical systems.

Other notes:

• Four regression tests are added.
• Only implemented in history-based mode.
• Currently does not support mesh-filtered tally with track-length estimator (due to the continuous weight capture). If such a tally is used, a warning shows up and automatically switches it to the collision estimator.
• Currently only use track length estimator for the alpha-related integral quantities. Collision and absorption estimators, similar to the keff tally, may be considered in future updates.
• When I manually used clang-format, I found that some files that I didn't touch got updated as well (e.g., lattice.*, string_utils.*, etc; I manually git-checked them out, though). Is that expected?

Thanks!