PP-LTG attempt 2 and fix related errors in adaptivity handling

Project.toml

@@ -1,6 +1,6 @@
 name = "OrdinaryDiffEqOperatorSplitting"
 uuid = "760fc936-9fa0-4281-9c1a-468957eedc89"
-version = "0.3.6"
+version = "0.4.0"
 authors = ["Dennis Ogiermann <termi-official@users.noreply.github.com> and contributors"]
 
 [deps]
@@ -19,7 +19,7 @@ Unrolled = "9602ed7d-8fef-5bc8-8597-8f21381861e8"
 [compat]
 BinaryHeaps = "1.0.0"
 CommonSolve = "0.2.4"
-DiffEqBase = "6.165.1, 7"
+DiffEqBase = "7.5"
 ModelingToolkit = "11"
 OrdinaryDiffEqCore = "4.4"
 OrdinaryDiffEqLowOrderRK = "1.7, 2"

_typos.toml

@@ -3,3 +3,5 @@
 Strang = "Strang"
 # Citation key for Trotter (Tro:1959:psg) used in docstrings
 Tro = "Tro"
+# Splitting variable name
+BA = "BA"

docs/Project.toml

@@ -6,4 +6,4 @@ OrdinaryDiffEqOperatorSplitting = "760fc936-9fa0-4281-9c1a-468957eedc89"
 [compat]
 Documenter = "1.16.1"
 DocumenterCitations = "1.4.1"
-OrdinaryDiffEqOperatorSplitting = "0.2.3, 0.3"
+OrdinaryDiffEqOperatorSplitting = "0.2.3, 0.3, 0.4"

docs/make.jl

@@ -31,6 +31,7 @@ makedocs(
         "Home" => "index.md",
         "usage/index.md",
         "Theory Manual" => "topics/time-integration.md",
+        "Adaptive time stepping" => "topics/adaptivity.md",
         "api-reference/index.md",
         "devdocs/index.md",
         "references.md",

docs/src/api-reference/index.md

@@ -17,6 +17,7 @@ GenericSplitFunction
 ```@docs
 LieTrotterGodunov
 StrangMarchuk
+PalindromicPairLieTrotterGodunov
 ```
 
 ## Per-node configuration

docs/src/topics/adaptivity.md

@@ -0,0 +1,96 @@
+# Adaptive time stepping
+
+Two independent layers of a splitting tree can adapt their step sizes:
+
+- **Leaf solvers** adapt their own internal steps within each sub-solve, exactly as
+  they would outside of this package. Their behavior is governed by the tolerances
+  and options they receive.
+- **Splitting nodes** adapt the splitting step itself. This requires an algorithm
+  that produces an error estimate of the splitting error;
+  [`PalindromicPairLieTrotterGodunov`](@ref) is currently the only one. It advances
+  with the average of the two mutually reversed Lie-Trotter sequences and uses half
+  their difference as the local error estimate.
+
+By default every node adapts exactly if its own algorithm can: for
+
+```julia
+alg = PalindromicPairLieTrotterGodunov((Tsit5(), Euler()))
+integrator = init(prob, alg; dt = 0.1)
+```
+
+the splitting node runs its step size controller, the `Tsit5` leaf adapts its inner
+steps, and the `Euler` leaf steps fixed. Passing `adaptive = false` (or a
+per-node [`TreeOption`](@ref)) overrides this.
+
+## The two tolerance layers
+
+`abstol` and `reltol` are understood at *every* node:
+
+- At a **splitting node** they scale the splitting error estimate: the node accepts
+  a step when `‖(u₁ - u₂)/2 ./ (abstol .+ max.(|u|) .* reltol)‖ ≤ 1`, where
+  `u₁, u₂` are the results of the two sequences of the palindromic pair.
+- At a **leaf** they control the accuracy of the inner solves.
+
+A scalar passed to `init` configures the whole tree with the same value. A
+[`TreeOption`](@ref) configures each node individually:
+
+```julia
+reltol = TreeOption(f, 1.0e-8)   # splitting node(s): tight
+reltol[1] = 1.0e-10              # leaf 1: tighter still
+reltol[2] = 1.0e-10              # leaf 2: tighter still
+integrator = init(prob, alg; dt, reltol)
+```
+
+## Rule: inner tolerances must be tighter than the splitting tolerances
+
+The palindromic error estimator measures the **difference of two inner-solved
+sequences**. Whatever error the inner solvers commit enters that difference as
+noise the estimator cannot distinguish from splitting error. Two consequences:
+
+1. **The estimator is blind to inner error.** With coarse fixed-step inner solvers
+   (say `Euler()` stepping at the splitting step size) the overall method degrades
+   to the inner order while the controller happily reports small `EEst`.
+2. **Loose inner tolerances choke the controller.** When a sub-problem has fast
+   internal dynamics (the typical reason for splitting it off), its solver commits
+   error on the order of its tolerance in *every* sub-solve, independent of the
+   splitting step size. If that exceeds the splitting tolerance — e.g. an inner
+   `abstol` 10× *larger* than the splitting `abstol` — then `EEst > 1` no matter
+   how small the splitting `dt` becomes. The controller keeps rejecting and
+   shrinking `dt` without the estimate improving, until the solve aborts with
+   `ReturnCode.DtLessThanMin` (or, with a tiny `dtmin`, grinds down to absurdly
+   small splitting steps).
+
+As a rule of thumb, keep the leaf tolerances at least **one to two orders of
+magnitude tighter** than the splitting tolerances (adaptive leaves), or the fixed
+inner steps well below the splitting step. The safe default is the scalar spread —
+identical tolerances everywhere are already borderline; never configure leaves
+*looser* than their splitting node.
+
+## Symptoms and causes
+
+| Symptom | Likely cause |
+|---|---|
+| `ReturnCode.DtLessThanMin`, `dt` collapsed, solution up to that point looks fine | Inner tolerances looser than the splitting tolerances (noise floor in the estimator), or genuinely unreachable tolerances |
+| Splitting `dt` grows to the full interval immediately | The operators (nearly) commute, the splitting error is ≈ 0; harmless — the inner solvers carry the accuracy |
+| Result visibly less accurate than the splitting tolerances suggest | Inner solves under-resolved: the estimator cannot see inner error (see rule above) |
+| A few rejections right after `init` or `reinit!` | Initial `dt` too large for the tolerance; harmless, the controller recovers |
+| Immediate abort with `DtLessThanMin` although tolerances look consistent | A scalar `dtmin` travels to the leaves too and may forbid the sub-steps they need. Restrict it to the splitting node with a `TreeOption` (`dtmin = TreeOption(f, 0.0); dtmin[] = 1e-3`) |
+
+## Controllers
+
+An adaptive splitting node runs an `OrdinaryDiffEqCore` step size controller
+(default: `IController`). The standard knobs (`qmin`, `qmax`, `gamma`,
+`qsteady_min`, `qsteady_max`, `failfactor`) can be passed to `init` — they are
+folded into the default controller — or a controller object can be passed
+explicitly via `controller` (per node via a `TreeOption`), e.g. a `PIController`
+whose memory smooths the step size sequence.
+
+## Failure handling
+
+A failing *adaptive* node (leaf or splitting node) is fatal: it already exhausted
+its own step size adaptation, and its return code propagates to the root. A failing
+*non-adaptive* node escalates the failure to the nearest adaptive ancestor, which
+rolls the whole subtree back and retries with a `failfactor`-shrunken step —
+shrinking the effective step of every non-adaptive descendant — until it either
+succeeds or falls below `dtmin`. Without any adaptive ancestor the integration
+stops with the escalated return code.

src/OrdinaryDiffEqOperatorSplitting.jl

@@ -15,7 +15,8 @@ import SymbolicIndexingInterface: variable_symbols
 import RecursiveArrayTools
 
 import OrdinaryDiffEqCore: OrdinaryDiffEqCore, isdtchangeable,
-    stepsize_controller!, step_accept_controller!, step_reject_controller!
+    stepsize_controller!, step_accept_controller!, step_reject_controller!,
+    accept_step_controller
 
 # In OrdinaryDiffEq v7 / DiffEqBase v7, passing verbose::Bool to inner ODE
 # integrators is no longer supported. Convert Bool → DEVerbosity when available.
@@ -50,7 +51,8 @@ include("integrator.jl")
 include("solver.jl")
 include("utils.jl")
 
-export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov, StrangMarchuk
+export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov, StrangMarchuk,
+    PalindromicPairLieTrotterGodunov
 export SplitNode, TreeOption
 
 include("precompilation.jl")

src/config_tree.jl

@@ -367,6 +367,30 @@ signed_dt_tree(config::ConfigTree, tdir, ::Type{tType}) where {tType} = ConfigTr
     map(child -> signed_dt_tree(child, tdir, tType), config.children)
 )
 
+"""
+    default_adaptive_option(f, alg)
+
+The default `adaptive` setting of `init`: a [`TreeOption`](@ref) in which every node
+adapts exactly if its own algorithm is adaptive. Passing a scalar `adaptive` instead
+configures the whole tree, including leaves whose algorithm cannot comply.
+"""
+function default_adaptive_option(
+        f::GenericSplitFunction, alg::AbstractOperatorSplittingAlgorithm
+    )
+    opt = TreeOption(f, SciMLBase.isadaptive(alg))
+    _default_adaptive!(opt, alg)
+    return opt
+end
+
+function _default_adaptive!(opt::TreeOption, alg)
+    opt.value = SciMLBase.isadaptive(alg)
+    inner = alg isa AbstractOperatorSplittingAlgorithm ? alg.inner_algs : ()
+    for (child, inner_alg) in zip(opt.children, inner)
+        _default_adaptive!(child, inner_alg)
+    end
+    return
+end
+
 """
     warn_non_adaptive(alg, config)
 
@@ -393,8 +417,14 @@ const NODE_OPTION_KEYS = (:dt, :adaptive, :verbose, :controller)
 inner_values(values::NamedTuple) =
     NamedTuple{filter(!in(NODE_OPTION_KEYS), keys(values))}(values)
 
-# ... of which a splitting node understands these.
-const SPLIT_OPTION_KEYS = (:dtmin, :dtmax, :failfactor, :isoutofdomain)
+# ... of which a splitting node understands these. The step-size controller knobs
+# (qmin, qmax, gamma, qsteady_min, qsteady_max) are not integrator options: they are
+# folded into the default controller of an adaptive node (`default_controller`) and
+# travel to the leaves like any other inner option.
+const SPLIT_OPTION_KEYS = (
+    :dtmin, :dtmax, :failfactor, :isoutofdomain,
+    :abstol, :reltol, :internalnorm,
+)
 
 function split_integrator_options(values::NamedTuple)
     inner = inner_values(values)