Allow the use of substeps for CRAM

[yrrepy]

Description

As part of my continuing battle against OpenMC deplete producing negative atom results and oscillations about zero.

This PR adds a substeps parameter to Integrator and SIIntegrator that subdivides each depletion interval into identical sub-intervals, reusing LU factorizations across them. This improves CRAM accuracy for nuclides with large λ·dt products — particularly those depleting toward zero, where single-step CRAM can undershoot to negative atom densities.

The implementation follows:
Isotalo & Pusa (2016), "Improving the Accuracy of the Chebyshev Rational Approximation Method Using Substeps," Nucl. Sci. Eng., 183:1, 65-77. https://doi.org/10.13182/NSE15-67

This addresses the same problem as #3879, but does it through improvements to solver numerical accuracy (step approximation error).
The two are complementary. Clip would be a safety net, and a nice cleanup of low atoms.

Substeps is invoked with:

openmc.deplete.PredictorIntegrator(operator, timesteps, power=power,
      substeps=2)  # subdivide each step into 2 sub-intervals

The default is left as substeps=1 (no substeps).

Personally, in future workflows with both PRs in, I would likely use both substeps=2, clip_min_atom_density=1e-20. As an improvement in numerical accuracy and a full-guard against negative atoms.

Similar script to #3879
demo.py

Here with substeps=2

[openmc.deplete] t=0.0 s, dt=1.0 s, source=100000000000000.0
[openmc.deplete] t=1.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=43201.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=86401.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=129601.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=172801.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=216001.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=259201.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=302401.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=345601.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=388801.0 s, dt=43200.0 s, source=0.0

WITH substeps=2 — improved accuracy near zero:
  step  1 (t=   0h):   5.79e+12
  step  2 (t=  12h):   1.78e-81
  step  3 (t=  24h):  5.47e-175
  step  4 (t=  36h):  1.68e-268
  step  5 (t=  48h):   0.00e+00
  step  6 (t=  60h):   0.00e+00
  step  7 (t=  72h):   0.00e+00
  step  8 (t=  84h):   0.00e+00
  step  9 (t=  96h):   0.00e+00
  step 10 (t= 108h):   0.00e+00
  step 11 (t= 120h):   0.00e+00

Attached in a zip (actually a .xz) are two further synthetic benchmarks that demonstrate the improvements.
substeps_synthetic-benchmarks.zip

They also demonstrate the near zero impact on calculation performance (with an MC inspired FOM).
The error is with respect to the case with fine explicit steps.
One can see from these results, despite improved behaviour, negatives that persist and could be cleaned up with the clip.

Checklist

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

[paulromano]

Thanks a lot for proposing this @yrrepy. I'm in favor of adding the ability to use substeps since it is well backed by the paper you mentioned and I see it is an option in Serpent. The PR conflates two changes right now though: adding substeps but also switching from spsolve to splu. This was originally proposed in #2703 and I think everyone was mostly on board. So I'd suggest we try to resolve any remaining issues with #2703 and get that merged first and then consider the substeps separately.

[yrrepy]

To riposte, the switch to splu only occurs for substeps>1, that is, within the new functionality; within a single solve. And splu here does not span step-to-step.

To frame it another way:
this PR implements intra-step splu (substeps), PR #2703 implements inter-step splu.
The inter-step splu cached LU factors across depletion timesteps with the same matrix. Where it encountered its hiccups.

Rather than a conflation, a different goal, substeps, which happens to benefit from a similar strategy.
The goals differ, no speed-up is desired or expected from the LU reuse here. This PR remains targeted and narrow in scope.

An alternative: proceed with this PR, and afterwards refactor the below to splu , with all the necessary inter-step piping

        if self.substeps == 1:
            A = dt * csc_array(A, dtype=np.float64)
            y = n0.copy()
            .....

[paulromano]

Understood, I did see that the use of splu only happens for substeps > 1. However, if we are going to use splu for substeps, I just want to make sure we understand whether there are any issues with splu in the first place, which would be resolved by handling #2703 first. I'm looking into that PR now so with any luck we should have both resolved soon 👍

[paulromano]

Quick update: see my comment on #2703 (review) regarding the use of splu. While it looks like it's not possible to use it across timesteps, there shouldn't be any issue with intrastep usage as you pointed out @yrrepy. I'll continue reviewing this one now!

[paulromano]

@yrrepy I did some refactoring here, mostly to change how substeps was passed through. To the user, it looks the same (as an argument on Integrator) but internally it is passed directly to the IFPCramSolver.__call__ method (rather than its constructor) to avoid having to create a separate instance for the sake of using different substeps. If you're good with the changes I made, I'd say we're good to merge!

[yrrepy]

Looks good to Perry, thanks for towing this across Paul!

It'd be great if we could get some others to test this and confirm for their cases that it reduces negatives (with and without the zero clip). And that performance impact is negligible.