3rd and 4th order methods - #109
Open
oscardssmith wants to merge 3 commits into
Open
Conversation
Implements the practical parts of Auzinger, Hofstaetter, Ketcheson & Koch, BIT 57:55-74 (2017): coefficient-table splitting schemes and the section 3 adjoint-pair local error estimator. Backward substeps. Every real splitting scheme of order three or above has negative coefficients, so some substeps run against the tree's direction -- previously unsupported, because a child's direction is fixed at construction and the stepping loop re-signs dt back before perform_step! sees it. reverse_direction! flips a leaf's tdir, dt, dtcache, dtpropose and dtmax and rebuilds its tdir-keyed heaps. A splitting node needs different handling: its heap stores raw times with the direction in the heap's type, so keys cannot be re-signed and only behind-times are dropped, and the reversal recurses into the subtree because add_tstop! propagates to descendants eagerly. Coefficient tables. SplittingCoefficients holds an s-by-N table as nested NTuples, checking each operator's consistency condition at construction. Ruth3 is the first table (exactly rational, order 3). One generic _perform_step! walks stages by operators and skips zero coefficients. LieTrotterGodunov and StrangMarchuk keep their hand-written steps. AdjointPair(base) averages a base scheme of odd order p with its adjoint for a solution of order p+1, and takes half their difference as the error estimate driving the controller. Even-order bases are rejected: their leading error terms are equal rather than opposite, so averaging would raise no order and the difference would not be an error estimate. The adjoint reverses the base's whole flat sequence of flows, reusing the same table. AdjointPair(LieTrotterGodunov(...)) is PalindromicPairLieTrotterGodunov, which a test asserts to the bit. Also fixes a pre-existing bug the new schemes exposed: the snap window absorbing child clock drift was scaled by the child's own dt, which after a substep can be arbitrarily small and collapse the window below the drift it must absorb. It now uses the outer step size, the scale the drift accumulated against. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Yoshida's triple jump, derived in this package's stage convention rather than transcribed: composing three Strang steps of lengths w1*h, w0*h and w1*h and merging the adjacent operator-1 flows gives four stages and eight flow evaluations rather than nine. The merge zeroes the last stage's second coefficient, so this is also the first table to exercise the zero-coefficient skip in the generic step. Measured order 4.02, 4.00, 4.00 over dt 0.2 to 0.025. w0 is negative and large: the second operator's cumulative time reaches 1.35*h before returning through -0.35*h, so this leans hard on backward substeps. Being of even order, Yoshida4 cannot serve as an AdjointPair base, which is now asserted alongside the StrangMarchuk case. Ruth3's two-operator restriction is factored into a shared helper, since Yoshida4's table is also AB only. Also adds the three bibliography entries the docstrings had been citing without defining -- Rut:1983:cim, Yos:1990:cho and AuzHofKetKoc:2017:psm -- which would have failed the docs build. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The citation key AuzHofKetKoc:2017:psm splits into Auz/Hof/Ket/Koc, and typos reads "Ket" as a misspelling of "Kept" or "Key". Added alongside the existing "Tro" exception for the Trotter key. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
termi-official
requested changes
Aug 4, 2026
termi-official
left a comment
termi-official
left a comment
Collaborator
There was a problem hiding this comment.
Looks already not bad! I left some comments. I think my main gripe is that we probably do not stress the rejection logic really hard right now. I think we should especially try to find ways to trigger step rejections mid way on some recursive split with different mixes of adaptivity and no adaptivity to see if the failure paths behave correctly and we do not get stuck in some faulty loops again. I had this problem in the past that a failing subintegrator was not handled correctly and we had to wail until the outer intergrators's dt starves slowly.
Comment on lines
+217
to
+223
| # This heap stores *raw* times and carries the direction in its ordering, which is | ||
| # part of its type, so unlike a leaf's its keys cannot be re-signed. | ||
| _drop_times_behind!(sub.tstops, sub.tdir, sub.t) | ||
|
|
||
| # `add_tstop!` propagates eagerly to every descendant, so a child left facing the | ||
| # old direction would reject the reversed node's next tstop as behind it. | ||
| _reverse_children!(sub.child_subintegrators) |
Collaborator
There was a problem hiding this comment.
I have trouble with this logic. Would it be easier to make the SplitSubIntegrator follow the convention in the other integrators?
| ) | ||
| parent.force_stepfail && return | ||
|
|
||
| # The average of the pair is the order p+1 solution ... |
Collaborator
There was a problem hiding this comment.
Should we reference to something in literature for this piece here, so future readers have a chance to understand where this is coming from?
Collaborator
There was a problem hiding this comment.
This file will grow quite quickly now when we start adding new schemes. Should we split it up and especially isolate the table logic?
Collaborator
There was a problem hiding this comment.
Devdocs on the tables is missing. Also, do we want to expand the theory section a bit for interested readers in the higher order theory and the connection with the tables?
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
based on https://link.springer.com/article/10.1007/s10543-016-0626-9