Fold RK46NL into the generic LowStorageRK2N dispatch - #3683
Conversation
|
family by family would probably be easier, yes. |
087bdb7 to
eedc675
Compare
|
Steps don't matter for non-adaptive methods, are the solutions identical? |
|
@ChrisRackauckas Yes, solutions are identical — validated by running each algorithm on a Lorenz problem and comparing the final state vector between the branch and master. |
| c5 = convert(T2, 0.8) | ||
| c2end = (c2, c3, c4, c5) | ||
|
|
||
| return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end) |
There was a problem hiding this comment.
these all already hit the same code, so it's not a substantive change. Do any other caches hit the same form?
|
Updated: helper indirection removed — Re 'do any other caches hit the same form': yes — |
|
@ChrisRackauckas i think once its merged we can do that then |
|
see the current comment. |
44cb756 to
fd240f0
Compare
|
@ChrisRackauckas ohh i see In previously merged |
fd240f0 to
5562000
Compare
singhharsh1708
changed the title
5562000 to
075cc7c
Compare
RK46NL was the last plain-2N method still carrying a bespoke constant cache struct, mutable cache struct, and four hand-unrolled initialize!/perform_step! methods. Its recurrence is the same Williamson form the shared LowStorageRK2N kernels already run, so the tableau moves into a RK46NLConstantCache function returning a LowStorageRK2NConstantCache with unchanged coefficients, RK46NL joins the LowStorageRK2NAlgorithm union, and the bespoke caches and kernels are deleted. The struct gains the williamson_condition keyword (default false) that the shared alg_cache reads, and isfsal/uses_uprev are declared false to match the rest of the family. The generic in-place kernel re-evaluates f at the start of each step instead of using FSAL storage, keeping the per-step f-call count at six. The new field is fourth, so positional construction now takes four arguments, the same break the nine sibling 2N methods already took when they gained the keyword. Fixed-dt end states are bit-identical to the removed implementation for in-place problems and for out-of-place problems whose state is an array. Out-of-place problems with a scalar state differ in the last bits: 3 ulp after 128 steps of u' = 1.01u, 1 ulp after 256 steps of a nonlinear scalar problem, roughly four thousand times below the method's own discretization error there. Both kernels expand to the same muladd(dt, k, A*tmp), but the straight-line code folds the other product into the fused multiply-add, and only scalar states are affected because muladd over arrays never fuses. RK46NL gains the williamson_condition testset the nine siblings already have. Version goes to 3.4.0 for the new public keyword.
075cc7c to
4a6bfc8
Compare
2 participants
Rework after master's LowStorageRK consolidation landed. The scaffold this PR originally proposed (a unified 2N tableau plus a new generic perform_step) was superseded by 570dde8, which introduced the LowStorageRK2NAlgorithm union, per-method tableau constructors, and shared alg_cache methods for nine 2N methods. RK46NL was the one plain-2N method left out: it still carried its own constant cache struct, mutable cache struct, and four hand-unrolled initialize!/perform_step! methods. This PR is now exactly that migration, following the same pattern as the earlier 2N and 3S consolidations.
The change is mechanical. The tableau moves into a RK46NLConstantCache function returning LowStorageRK2NConstantCache with unchanged coefficients, RK46NL joins the LowStorageRK2NAlgorithm union, and the bespoke caches and kernels are deleted. The struct gains the williamson_condition keyword (default false) that the shared alg_cache reads, and isfsal/uses_uprev are declared false to match the rest of the family. The generic in-place kernel re-evaluates f at the start of each step instead of using FSAL storage, so the per-step f-call count stays at six. Net -134 lines of source.
One thing to be aware of: williamson_condition is now the struct's fourth field, so RK46NL(stage_limiter!, step_limiter!, thread) no longer resolves and positional construction takes four arguments. That is the same tradeoff the nine sibling 2N methods already made when they gained the keyword, and keyword construction is unchanged.
Version goes to 3.4.0. Master released OrdinaryDiffEqLowStorageRK 3.3.0 in #3241, and this adds a public keyword argument to an exported algorithm struct.
The previous revision of this branch was numerically broken in place: it declared isfsal(RK46NL) = false while keeping a kernel whose first stage read cache.fsalfirst, which is only refreshed by the FSAL copy, so every step after the first used a stale derivative and measured order collapsed to ~0. Deleting that kernel in favor of the shared one removes the failure mode instead of patching it.
A correction to the earlier text of this description, which claimed fixed-dt end states are bit-identical to the removed implementation. That holds for most but not all of the space. Running one script against master at 0be95df and against this branch, dt = 1/128 and adaptive = false: every in-place problem tested comes out bit-identical (2x2 linear, Lorenz, a nonlinear scalar wrapped in a one-element vector), and so does every out-of-place problem whose state is an array (2x2 linear, out-of-place Lorenz). Out-of-place problems whose state is a plain Number do not. Those land 3 ulp apart after 128 steps of u' = 1.01u (4.9e-16 relative) and 1 ulp apart after 256 steps of u' = 1.01u - 0.3u^2 + 0.1sin(t) (1.7e-16 relative). On the linear problem both trees sit 1.9e-12 from the analytic solution, so the disagreement between them is about 3900 times smaller than the method's own discretization error at that step size.
The cause is fused multiply-add contraction, which is also why only scalars are touched. Both the removed kernel and the shared one write the stage update as Atmp + dtk, and @muladd expands both to the same muladd(dt, k, Atmp). Compiled, the straight-line unrolled version folds the other product into the fma instead: putting fma(A, tmp, dtk) into a standalone replica of the recurrence reproduces master's bits exactly on both scalar problems, and fma(dt, k, A*tmp) reproduces this branch's. When the state is an array, muladd is not a fused operation at all, both forms round the same way, and every array case matches to the bit.
Testing. RK46NL had no williamson_condition coverage, so this adds an RK46NL testset next to ORK256 built the way the nine siblings are: williamson_condition = true and false both go through test_convergence on the only-time, linear, and nonlinear problem sets at dts = 1 ./ 2 .^ (7:-1:3), followed by the Base.summarysize register bounds and the alias_u0 check. Measured order is 3.97 on the only-time problems, 3.99 on the three linear ones and 4.03 on the nonlinear ones, identical for both settings of the keyword. The full ode_low_storage_rk_tests.jl file now reports 566 pass, 28 pre-existing broken and 0 failures, against 547 pass and the same 28 broken on master, so the 19 new tests all pass and nothing else moved. The 515/28 figure quoted earlier in this description predates the rebase onto #3241. AllocCheck finds perform_step! allocation-free for RK46NL on the new cache path with williamson_condition either way.
AI Disclosure
Claude assisted with this work.