openmc-dev
Adding explicit LU factorization within CRAM
Description
Another performance (and accuracy?) improvement to depletion. For some reason doing an explicit LU factorization before solving within the CRAM method seems to result in far fewer warnings about negative density nuclides and gives a ~10% speed up in the cases I have tried. The testing isn't exhaustive so I'd be curious to find out what others observe. I have attached files below showing the output during the depletion calculation and profile information for my particular example problem. I don't yet fully understand these differences since it gets into the scipy and SuperLU internals and may ultimately end up being very machine dependent because of this, but I'll keep digging and if anybody has any insight, feel free to chime in.
Update after merging #2764: no_patch.txt with_patch.txt
There is almost a 2x speed up for problems with many repeated timesteps with the same source rate
Checklist
- [x] I have performed a self-review of my own code ~~- [ ] I have run clang-format on any C++ source files (if applicable)~~
- [x] I have followed the style guidelines for Python source files (if applicable)
- [x] 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)~~
Thanks for taking a look @drewejohnson! All good points.
Now, the question becomes does this change make the results more correct, and thus the failure is because the reference files need to be updated?
I have tried updating the reference results and they still don't pass in CI, which I have struggled to figure out why and wonder if it has to do with how scipy is built and/or linked on the CI servers vs my own machine.
It would be great if you could add a test that shows negative densities w/ the current implementation, and fewer or no negative densites with superlu.
Anecdotally you can see the results in the output of the two text files I posted. The no_patch version has a number of negative density warnings and the with_patch version does not. The code as written doesn't give the user the ability to run depletion calculations with the current implementation and the new implementation so a formal test case isn't possible (if that's what you were suggesting).
If I'm a user working on a depletion model, and there are two linear solver methods for depletion, how do I know what to use and why?
The user never really interacts with the CRAM method specifically and I would like to get rid of this kwarg anyway, but I can't figure out how to get these tests to pass both locally and in CI without this logic to revert to the direct solve when transfer_rates is set on the Integrator.
FWIW the tests in tests/regression_tests/deplete_with_transfer_rates/ do pass locally for me on this branch since adding the direct_solve kwarg to keep the behavior the same as before when transfer_rates is set.
I have updated the profile results for this patch (difference between current develop and this branch) since some of the other improvements have been made. This is the same depletion problem as original posted at the top of this PR. There is still a ~10% speed up and fewer negative density warnings as expected.
I don't love the extra logic in the CRAM method here, but I think this is the cleanest approach for now unless we would prefer to create a separate solver class for use when transfer_rates is set. @paulromano and @drewejohnson what do you think?
I'm in favor of the LU factorization because of the alleged performance boost. Feels great.
it's not immediately obvious that the assumptions made in the CRAM method are satisfied when trying to solve the block system that arises from the use of transfer rates
Yeah this feels unrelated to using csc / csr / lu factorization thing at hand: Is CRAM appropriate w/ transfer rates. I wasn't around for the transfer rates discussion so I'm going to assume its related to moving fuel or depleting systems where an external actor is supplying / removing nuclides. I keep going back to this section from Maria Pusa's paper on CRAM
However, it was recently discovered by the author that the eigenvalues of the burnup matrix are generally confined to a region near the negative real axis.1 This observation led to applying the Chebyshev Rational Approximation Method to solve the burnup equations. This method can be interpreted as the best rational approximation on the negative real axis
So if you're using transfer rates as a way to scrub fuel, like off gassing fission products, then the diagonal of the matrix will remain negative. Probably okay. But, if you include enough addition of material such that the you have a positive rate of change, CRAM may not be valid.
I personally don't like the inclusion of the conditionals in the CRAM solve. Is the rationale behind not using LU for transfer rate problems we can't get them to pass on CI? Is there a way to generate the reference data on the CI machine and pull those? Otherwise we're at the mercy of machine differences, library versions, compilers, all those fun things python hides from us :slightly_smiling_face:
I personally don't like the inclusion of the conditionals in the CRAM solve. Is the rationale behind not using LU for transfer rate problems we can't get them to pass on CI? Is there a way to generate the reference data on the CI machine and pull those?
See PR #2764 @drewejohnson I think I'm going to table this one until that one is sorted out. I think once we understand the issues in the tests with that PR then this one will sort itself out. In that PR I didn't change the solver, just the matrix format and the tests still didn't pass. I had to loosen the tolerance significantly for them to pass. I would like to still find out the fundamental cause, but I think it's unrelated to the splu object as evidenced in #2764.