Methods and step-size control

Status: M1 draft, 2026-10-23. Normative for both implementations (D1, D8) unless a paragraph is marked informative or assumption.

This chapter defines what one step of each v1 method computes in the forward pass, how its error is estimated and measured, and how a leaf chooses, accepts and rejects its steps. The discrete adjoint of the same steps is adjoint.md ; the two chapters use the same arithmetic forms, and this one refers there for the step equations rather than repeating them. What the run config may set is run-config.md §4.2.

1. Sources and the coefficient tables

What is borrowed (D31). Coefficients, error estimators, controller formulas and controller constants are OrdinaryDiffEq’s, read from these installed versions:

KPP’s Rosenbrock methods come from KPP 2.1 as shipped in WRF-Chem (WRF/chem/KPP/kpp/kpp-2.1/int/WRF_conform/rosenbrock.f90, relative to the directory holding the repositories). WRF_RK3 comes from Wicker and Skamarock (2002) and WRF’s dyn_em/solve_em.F:610-622. IMEXEuler and Theta are defined here and in adjoint.md §4.4–§4.5.

The tables are data. Every method’s coefficients, error estimator, order, flags and controller constants live in tableaus/<name>.json , one file per run-config alg value, validated by schemas/tableau.schema.json . Both implementations MUST read their coefficients from these files, or generate code from them; neither restates the numbers by hand. Each file cites its source by package, version, file and line.

The files are generated by ../tools/tableaus/extract.jl , which reads OrdinaryDiffEq’s tableau constructors where they exist and transcribes, with a citation, the coefficients OrdinaryDiffEq writes inline in perform_step! (Heun, RK4, ForwardEuler and the three SSP methods). Before it writes, the script checks every table: one step evaluated from the table in the form this chapter defines matches OrdinaryDiffEq’s own step to roundoff (to the Newton tolerance for SDIRK), the error estimate matches OrdinaryDiffEq’s where it is a plain embedded estimate, and every method shows its order on a smooth non-autonomous test problem. Regenerating produces byte-identical files.

Two tables deliberately differ from OrdinaryDiffEq (§5.4, §5.5); their files say so in notes.

2. Method catalogue

“Error” is the error_estimate field of the table. Methods with error none have no estimate, so a leaf using them MUST give dt (D52). Orders are OrdinaryDiffEq’s alg_order and alg_adaptive_order; KPP’s are from KPP’s ros_ELO.

Informative. ROS2 and KPP_Ros2 show an observed order of 1.75 at the step sizes the extraction check uses, and Kvaerno4 3.70, which is pre-asymptotic. The others are within 0.15 of their order.

3. Explicit Runge–Kutta

3.1 Forms and evaluation order

Butcher form (A, b, c): the equations of adjoint.md §4.1. Stage $i$ is evaluated at time $t_n + c_i h$. Each stage argument and $u_{n+1}$ is summed left to right starting from $u_n$, then over ascending $j$, each term formed as $(h,a_{ij}),K_j$ with the scalar product first. These sums are not contracted into fused multiply-adds (D10). Deviation: OrdinaryDiffEq writes its steps under @muladd, which lets the compiler fuse, and groups the stage sums as $u_n + h,(\sum_j a_{ij} K_j)$, so it agrees with this form to roundoff, not bit for bit. A zero coefficient contributes nothing (the term is skipped, not multiplied by zero), so a non-finite $K_j$ never reaches a stage that does not read it.

Shu–Osher form (alpha, beta, c), for SSPRK33, SSPRK43 and SSPRK104. With $u^{(0)} = u_n$,

u^{(i)} = \sum_{j=0}^{i-1} \Big(\alpha_{ij}\, u^{(j)} + (h\,\beta_{ij})\, f\big(u^{(j)}, t_n + c_j h\big)\Big), \qquad u_{n+1} = u^{(s)} ,

summed flat over ascending $j$, the $\alpha$ term before the $\beta$ term for each $j$, skipping zero coefficients. alpha[i-1][j] and beta[i-1][j] hold $\alpha_{ij}$ and $\beta_{ij}$, and c[j] is the time of the evaluation at $u^{(j)}$; each $f(u^{(j)})$ is evaluated once. OrdinaryDiffEqSSPRK writes these steps inline in a regrouped arithmetic (for SSPRK104, (3u0 + 2u4 + 2(h/6)k4)/5), which agrees with this form to roundoff but not bit for bit; the tables are the regrouped coefficients (verified to $\le 2\cdot10^{-16}$ relative by the extraction check).

Interpolation and dense output are not part of v1 (output times are tstops, run-config §7.3), so the tables carry no interpolation coefficients.

3.2 Error estimates

• embedded (Heun, BS3, Tsit5): $e = \sum_i (h,\tilde b_i), K_i$ with $\tilde b$ = btilde, the scalar product first. Deviation: OrdinaryDiffEq forms $h,(\sum_i \tilde b_i K_i)$ (OrdinaryDiffEqTsit5 src/tsit_perform_step.jl:167). In BS3 and Tsit5 the last stage is $f(u_{n+1}, t_n + h)$ (row $s$ of $A$ equals $b$), so the estimate reads it.

• shu_osher_embedded (SSPRK43): the embedded solution $\tilde u = \sum_{j=0}^{s} \tilde\alpha_j u^{(j)} + h\sum_{j=0}^{s-1}\tilde\beta_j f(u^{(j)})$ with embedded.alpha and embedded.beta, and $e = \sigma(\tilde u - u_{n+1})$ with $\sigma$ = embedded.scale. For SSPRK43, $\tilde u = \tfrac13 u_n + \tfrac23 u^{(2)} + \tfrac13 h f(u^{(2)})$ and $\sigma = \tfrac12$ (OrdinaryDiffEqSSPRK src/ssprk_perform_step.jl:766-807).

• residual_control (RK4), Shampine’s residual estimate (OrdinaryDiffEqLowOrderRK src/fixed_timestep_perform_step.jl:252-316). With $k_1 = f(u_n, t_n)$, $k_5 = f(u_{n+1}, t_n + h)$, and for each $\sigma \in$ residual_control.sigma $= \tfrac12 \mp \tfrac{\sqrt 3}{6}$:

  p_\sigma = (1-\sigma)u_n + \sigma u_{n+1} + \sigma(\sigma-1)\big[(1-2\sigma)(u_{n+1}-u_n) + (\sigma-1)h k_1 + \sigma h k_5\big],
  

  p'_\sigma = k_1 + \frac{\sigma}{h}\big[-4hk_1 - 2hk_5 - 6u_n + \sigma(3hk_1 + 3hk_5 + 6u_n - 6u_{n+1}) + 6u_{n+1}\big],
  

  $E_\sigma = \lVert h,(f(p_\sigma, t_n + \sigma h) - p’\sigma)\rVert$ in the norm of §7, and $E = 2.1342,\max(E{\sigma_1}, E_{\sigma_2})$ (residual_control.scale). This costs two extra $f$ evaluations per step. Unlike the other estimates it is a scaled norm already, so §7’s norm is applied to each residual, not to a single error vector.

• none: no estimate; the leaf runs at fixed dt.

3.3 First same as last

For the methods whose table has first_same_as_last: true (Tsit5, BS3, Heun, RK4, ForwardEuler, Rosenbrock23), the step’s last evaluation is $f(u_{n+1}, t_n + h)$, and the next step’s first stage is that value. It is evaluated at $t_n + h$ even when the step’s end time was set to a stop exactly (§10), as OrdinaryDiffEq does; replay and the adjoint use the time it was actually evaluated at. A leaf MAY reuse it only while $u$, $t$, the parameters and the resident forcing are unchanged since it was computed; setting $t_{n+1}$ to the stop exactly does not count as a change of $t$. A parameter bound by bind_dt: "leaf" changes with the step size, so on a leaf with one the value is reused only by a step whose size equals, bit for bit, the size of the step that computed it. It MUST evaluate afresh:

• at entry to every invocation of the leaf, because maps and other leaves may have changed the state (adjoint.md §4.1);
• after a domain clip (§11), after a skipped tiny stop step (§10), and after a forcing epoch change.

After a rejected step the start point is unchanged, so the stored first stage stays valid. ForwardEuler needs the end-of-step evaluation only for this reuse: a leaf MAY skip it and evaluate $f(u_n)$ at the start of each step instead.

4. Rosenbrock

4.1 Step

• Transformed form (A, C, gamma, c, d, b; ROS2, ROS3P, Rodas3/4/5P and the KPP methods): the equations of adjoint.md §4.2, $\hat W = I/(h\gamma) - J(u_n)$, stage $i$ evaluated at $t_n + c_i h$ with $h,d_i,f_t(u_n)$ added. This is OrdinaryDiffEqRosenbrock’s generic step (src/rosenbrock_perform_step.jl:431-540), which stores $W = J - I/(h\gamma)$ and negates the solve.
  • Arithmetic. $\hat W$’s diagonal entries are $1/(h\gamma) - J_{ii}$, each rounded once; its off-diagonal entries are $-J_{ij}$. The stage right-hand side is grouped $r_i = \big(f(U_i) + (h d_i), f_t\big) + s_i$, with $s_i = \sum_j (C_{ij}/h), k_j$ summed in ascending $j$ starting from zero; the $f_t$ term is absent when $d_i = 0$ or $f_t$ is not formed (§4.2). This is the grouping of OrdinaryDiffEqRosenbrock 2.7.5’s out-of-place step (:487-497); its in-place _stage_rhs! (:645-663) groups differently. $U_i = u_n + \sum_j A_{ij} k_j$, $u_{n+1} = u_n + \sum_i b_i k_i$ and $e = \sum_i \tilde b_i k_i$ are summed in ascending order. Every sum here skips zero coefficients, where OrdinaryDiffEq multiplies them through; the two differ at most in the sign of a zero.
• Rosenbrock23 (d, c32): the ode23s step of adjoint.md §4.2 (OrdinaryDiffEqRosenbrock src/rosenbrock_perform_step.jl:249-330). $W = I - h d,J(u_n)$ mathematically; in arithmetic, $M = I/(h d) - J(u_n)$ is factored and $W^{-1} r$ is computed as $(M^{-1} r) \cdot \tfrac{1}{h d}$, as OrdinaryDiffEq does with the sign moved. The adjoint’s transposed solves are scaled the same way. With $\sigma = 1/(hd)$ the step is $k_1 = W^{-1}(f_0 + h d, f_t)$, $U_2 = u_n + \tfrac{h}{2} k_1$ at $t_n + h/2$, $k_2 = \big(M^{-1}(f_1 - k_1)\big)\sigma + k_1$, $u_{n+1} = u_n + h k_2$. $f_2 = f(u_{n+1}, t_n + h)$ is evaluated on every step, fixed steps included, because it is the next step’s $f_0$ (§3.3); $k_3$ needs no further evaluation, so a leaf makes two right-hand-side calls per step after its first, adaptive or not (41 for forced_decay at 20 fixed steps).
• From the third stage on, a stage whose row of $A$ and whose $c_i$ equal the previous stage’s MAY reuse the previous stage’s $f$ value. OrdinaryDiffEq does this (:463-466), which applies to ROS3P, KPP_Ros3 and KPP_Ros4, not to Rodas3 or Rodas4; KPP flags it in kpp.new_f. This changes no arithmetic, only right-hand-side counts.

4.2 Jacobian, time gradient and $W$

• $J$ and $f_t$ are formed at $(u_n, t_n)$ on every step attempt, from the leaf’s band models (the jacobian option), and $W$ is factored once per attempt. v1 Rosenbrock leaves never reuse $J$ across steps (adjoint.md §4.2).
• Informative. OrdinaryDiffEq reuses $J$ across steps for W-methods (OrdinaryDiffEqDifferentiation src/derivative_utils.jl:46-171). Rosenbrock23 is marked a W-method but defaults to max_jac_age = 1, which forms $J$ every step; the Rodas and ROS methods are not W-methods and always form $J$. So the forward behaviour matches OrdinaryDiffEq’s defaults for every v1 Rosenbrock method.
• $f_t$ is formed only when some $d_i \ne 0$ and the leaf’s operators depend on $t$ explicitly or through forcing within the current epoch. Otherwise it is zero.
• A singular $W$ (a failed factorization) always rejects the step, and the new size is the controller’s rejection size at $E = 2$, as OrdinaryDiffEq sets set_EEst!(integrator, 2) (src/rosenbrock_perform_step.jl:272, :448) (§8.5). The rejection holds under the PID controller too, whose factor at $E = 2$ can exceed pid_accept_safety. At fixed dt a singular $W$ ends the leaf with Unstable (§12).
• Linear solves. dense_lu and block_lu factor by Julia’s generic_lufact! arithmetic: partial pivoting choosing the first row of largest magnitude (strict comparison), multipliers $a_{ik} \cdot (1/a_{kk})$, and column-oriented triangular solves (forward over ascending columns, back over descending columns with $x_j = y_j / u_{jj}$). The transposed solve of the adjoint uses the same factors. block_lu applies this to every lane at once; with pivot: false it skips the row exchanges. A LAPACK factorization would differ only at roundoff, but the two implementations agree bitwise only under this rule.
• Precision. $W$ and its factors may be stored in a lower precision than the state (architecture §8, the jacobian and factors roles). Classical Rosenbrock methods assume an exact $J$, so the planner refuses a low-precision $W$ with them unless overridden (run-config §7.1).

4.3 Error estimates

• Transformed form, embedded: $e = \sum_i \tilde b_i k_i$ with $\tilde b$ = btilde (no factor $h$: the $k_i$ already carry it).

• Rosenbrock23: with $f_0 = f(u_n, t_n)$, $f_1 = f(U_2, t_n + h/2)$ and $f_2 = f(u_{n+1}, t_n + h)$,

  W k_3 = f_2 - c_{32}(k_2 - f_1) - 2(k_1 - f_0) + h d\, f_t(u_n), \qquad e = \tfrac{h}{6}\,(k_1 - 2k_2 + k_3) .
  

  $f_2$ is the next step’s $f_0$ (§3.3).

  The $f_t$ term is $h d, f_t$, the same as the first stage’s, which makes $e$ a third-order estimate, $O(h^3)$, on problems that depend on $t$ explicitly (D62). Deviation: OrdinaryDiffEq adds $h f_t$ (src/rosenbrock_perform_step.jl:104, :298). That leaves $\tfrac{h}{6}, h (1 - d), f_t$ in the estimate, which is then second order (1.99 on forced_decay, against 3) and overstates the error as $h \to 0$, so its steps on time-dependent problems are smaller than the tolerance requires. $u_{n+1}$ is the same either way, so fixed-step runs and the oracle (D31) are unaffected; adaptive Rosenbrock23 step sequences on time-dependent problems differ from OrdinaryDiffEq’s.

4.4 KPP methods

KPP’s ros_A, ros_C, ros_Gamma(1), ros_Alpha, ros_Gamma, ros_M, ros_E are this form’s $A$, $C$, $\gamma$, $c$, $d$, $b$, $\tilde b$: KPP solves $(I/(h\gamma) - J)K_i = f(Y_i) + \sum_j (C_{ij}/h) K_j + h\gamma_i f_t$ (rosenbrock.f90:575-627). KPP’s error norm is §7’s. Its step-size rule differs: an I controller with exponent $1/$ros_ELO, safety 0.9, bounds $[0.2, 6]$, no growth after a rejection, and a factor 0.1 after two rejections in a row (rosenbrock.f90:631-660); the table records it under kpp.controller. By default a KPP leaf uses §8’s controller with the constants in its controller field.

KPP’s own step-size rule becomes a run-config controller option, controller: "KPP", at M7 for WRF-Chem parity (D58). Until then KPP methods run under the PI controller.

5. SDIRK and ESDIRK

5.1 Stage equations

The tables give $A$, $b$, $c$ with $a_{ii} = \gamma$ on implicit stages. The unknown of stage $i$ is $z_i = h K_i$:

\tau_i = u_n + \sum_{j<i} a_{ij}\, z_j, \qquad
z_i = h\, f\big(\tau_i + \gamma z_i,\ t_n + c_i h\big), \qquad
u_{n+1} = u_n + \sum_i b_i\, z_i ,

which is adjoint.md §4.3 with $Z_i = \tau_i + \gamma z_i$. This is OrdinaryDiffEqSDIRK’s generic ESDIRK step (src/generic_imex_perform_step.jl:1400 onward, the out-of-place body).

Arithmetic (D10). $\tau_i$ and $u_{n+1}$ are summed left to right starting from $u_n$, over ascending $j$ (or $i$), each term $a_{ij}, z_j$ (or $b_i, z_i$), skipping zero coefficients, without fused multiply-adds. $u_{n+1}$ is always $u_n + \sum_i b_i z_i$, never $\tau_s + \gamma z_s$, also for the stiffly accurate methods, where OrdinaryDiffEq takes the latter; the two agree to roundoff.

An explicit first stage (explicit_first_stage: true: TRBDF2, Kvaerno3, Kvaerno4) is $z_1 = h,f(u_n, t_n)$, evaluated, componentwise $h \cdot f_j$. OrdinaryDiffEq instead reuses the previous step’s $z_s/h$ (src/generic_imex_perform_step.jl:2422-2448), which equals $f(u_n)$ only to the Newton tolerance. The adjoint (adjoint.md §4.3) differentiates $K_1 = f(u_n)$, so the forward pass must compute exactly that.

5.2 Newton iteration

Each implicit stage is solved by modified Newton with $W = I - h\gamma J$ (§5.3), following OrdinaryDiffEqNonlinearSolve’s NLNewton (src/nlsolve.jl:105-275, src/newton.jl:708-756, defaults src/type.jl:185-190):

1. Initial guess. $z_i^{(0)} = z_{i-1}$ for implicit stages after the first, and $z^{(0)} = 0$ for the first implicit stage of a step and for the single implicit stage of §6. Assumption: OrdinaryDiffEq’s per-method predictors (Predictor.Linear and stage extrapolation) change only the iteration count, not the converged stage beyond the Newton tolerance, so v1 specifies this simpler guess.

2. Iteration $m = 1, 2, \dots$: $\ \Delta z = W^{-1}\big(h f(\tau_i + \gamma z^{(m-1)}, t_n + c_i h) - z^{(m-1)}\big)$, $\ z^{(m)} = z^{(m-1)} + \Delta z$.

   • Arithmetic. With $\sigma = 1/(h\gamma)$ and $M = \sigma I - J$ (§5.3), componentwise: $Z_j = \tau_{i,j} + \gamma z_j$, $r_j = (h, f_j(Z, t_n + c_i h) - z_j)\cdot\sigma$; then $\Delta z = M^{-1} r$ and $z_j \leftarrow z_j + \Delta z_j$. This is OrdinaryDiffEqNonlinearSolve’s _compute_rhs and dz = W \ ztmp (src/newton.jl:874-893, :735) with its $W = J - \sigma I$ negated, so $z^{(m)}$ equals its z .- dz bit for bit.

3. Increment norm $\delta_m = \lVert \Delta z \rVert$ in §7’s norm, scaled with $\max(|u_n|, |\tau_i + \gamma z^{(m-1)}|)$. A fixed-dt implicit leaf uses its control tolerances for this norm too (run-config §4.2).

4. Rate. For $m > 1$, $\theta_m = \max(0.3,\theta_{\text{prev}}, \delta_m/\delta_{m-1})$; for $m = 1$, $\theta_1 = \theta_{\text{prev}}$. $\theta_{\text{prev}}$ is 1 at the start of every step attempt and carries across that attempt’s stages only. It is set to $\theta_m$ when iteration $m$’s $\Delta z$ is applied, not when step 6 or 7 stops the iteration first (nlsolve! stores prev_θ after those tests, src/nlsolve.jl:228-230). Deviation: OrdinaryDiffEq carries $\theta_{\text{prev}}$ across steps and resets it to 1 only in the attempts that follow a rejection (:222-224). The tape does not record where rejections fell (spec/tape.md §2), so replay could not reproduce such a step bitwise (D10); here the step depends only on $(u_n, h)$, as §5.3 intends (D63).

5. Convergence when $m = 1$ and $\delta_1 < 10^{-5}$, or when $\eta_m \ge 0$ and $\eta_m,\delta_m < \kappa$, with $\eta_m = \theta_m/(1 - \theta_m)$ on every iteration, the first included, and $\kappa = 1/100$. nlsolve! sets an $\eta$ from $\eta_{\text{old}}$ before the loop (:140) but overwrites it on iteration 1 before any test (:235), so $\eta_{\text{old}}$ has no effect on convergence and is not kept. With $\theta_1 = 1$, $\eta_1$ is infinite and the $\eta$ test fails; with $\theta_1 > 1$, $\eta_1$ is negative and it fails. On a fixed-dt leaf the $\eta$ test is skipped on iteration 1 (check_η_convergence, :238-251), so there iteration 1 converges only through $\delta_1 < 10^{-5}$.

6. Floating-point limit. For $m > 1$, if $|\theta_m - 1| \le 100\sqrt\epsilon$ (eps_around_one, $\epsilon$ of f64): converged when $\delta_m \le 1$, diverged otherwise.

7. Divergence when $m > 1$ and $\theta_m > 2$ (strictly), or when $\delta_m$ is not finite.

   Steps 6 and 7 are tested before step 5. When either stops the iteration, the iteration’s $\Delta z$ is not applied, as in nlsolve!: a stage converged at the floating-point limit is $z^{(m-1)}$.

8. Failure after 10 iterations without convergence, or on divergence. The step is rejected and $h \leftarrow h/2$ (failfactor, OrdinaryDiffEqCore src/integrators/controllers.jl:488-491). A fixed-dt leaf cannot shrink its step, so a Newton failure there ends the leaf with ConvergenceFailure (§12).

The same iteration solves the implicit part of BackwardEuler, Theta and IMEXEuler (§6), and the converged stages are what the adjoint differentiates (adjoint.md §4.3, D42).

5.3 Jacobian and $W$

$J$ is formed at $(u_n, t_n)$ on every step attempt, and $W = I - h\gamma J$ is factored once per attempt and shared by all implicit stages (SDIRK has one $\gamma$). In arithmetic the factored matrix is $M = \sigma I - J(u_n, t_n)$ with $\sigma = 1/(h\gamma)$, $\sigma$ rounded once from the rounded product $h\gamma$, and $M$ formed as §4.1 forms $\hat W$, and $W^{-1} r$ is computed as $M^{-1}(r \cdot \sigma)$ (§5.2). Linear solves, a singular $M$ (a rejection, or Unstable at fixed dt) and precision roles are as in §4.2.

Informative. OrdinaryDiffEq reuses $J$ across steps and rebuilds $W$ only when $h\gamma$ changes by more than 20 % (new_W_dt_cutoff = 1/5), after a divergence with a stale $J$, or after an error-test failure on the first stage (OrdinaryDiffEqDifferentiation src/derivative_utils.jl:541-600). v1 does not: forming $J$ every attempt makes the step a function of $(u_n, h)$ alone, which keeps lane pools (D17) and replay (architecture §10) simple. The adjoint is unaffected either way (D42). Adopting the reuse rule for large column solves is revisited at M4b (D57).

5.4 Stage times

Stage $i$ is evaluated at $t_n + c_i h$ with $c_i = \sum_j a_{ij}$, the row sum the method’s order conditions assume. The tables store these.

OrdinaryDiffEqSDIRK sets $c_2 = \gamma$ for Kvaerno3 and Kvaerno4, where the row sum is $2\gamma$ (src/imex_tableaus.jl:1203-1272, :1275-1360, commented “matching master’s nlsolver.c = γ”). The tables keep OrdinaryDiffEq’s value as c_ordinarydiffeq. The extraction check shows that OrdinaryDiffEq does evaluate stage 2 at $\gamma$, and on the non-autonomous test problem Kvaerno3 and Kvaerno4 then converge at first order (observed 0.99 and 0.98) instead of third and fourth. With the row sums they converge at 2.94 and 3.85. Autonomous problems are unaffected.

Tier-2 trajectory comparisons against OrdinaryDiffEq’s Kvaerno3/4 use autonomous problems only (D55).

5.5 Error estimates

• embedded_smoothed (TRBDF2, SDIRK2, Kvaerno3, Kvaerno4): $e_{\text{raw}} = \sum_i \tilde b_i z_i$, smoothed as Hosea and Shampine (1996) propose:

  e = (I - h\gamma J)^{-1}\, e_{\text{raw}} ,
  

  one extra solve with the step’s $W$, computed as $e = M^{-1}(e_{\text{raw}} \cdot \sigma)$ (§5.3). $e_{\text{raw}}$ is summed over ascending $i$, each term $\tilde b_i, z_i$, skipping zero coefficients.

  Deviation from OrdinaryDiffEq. OrdinaryDiffEqSDIRK smooths with its stored Newton matrix $W_{\text{ODE}} = J - I/(h\gamma)$ (src/generic_imex_perform_step.jl:2389-2392), which gives $-h\gamma,(I - h\gamma J)^{-1} e_{\text{raw}}$: the estimate shrinks by the factor $h\gamma$. A run on $f(u, t) = t^3$ (where $J = 0$ and smoothing should change nothing) confirms it: OrdinaryDiffEq’s smoothed TRBDF2 estimate is $h\gamma$ times its unsmoothed one (ratio 34.1 at $h = 0.1$, 341 at $h = 0.01$, with $\gamma = 0.2929$). With smooth_est = false OrdinaryDiffEq’s estimate equals $e_{\text{raw}}$, which the extraction check confirms. The spec adopts the formula above (D56). The rule is to follow OrdinaryDiffEq unless that gives wrong or low-accuracy results, and here it does: the estimate’s extra factor $h\gamma$ makes the controller take steps whose error grows far beyond the tolerance as the tolerance tightens. ../tools/tableaus/sdirk_error_estimate_check.jl measures the end-of-run error divided by the tolerance (reltol = abstol) for OrdinaryDiffEq with and without smoothing, and for this section’s step (with the row-sum stage times of §5.4) under §8.2’s PI controller with the estimate above:

  Problem	Method	tol	smoothed (OrdinaryDiffEq default)	unsmoothed	adopted
  forced_decay	TRBDF2	1e-6	165	9.4	8.7
  forced_decay	TRBDF2	1e-8	1430	40	43
  forced_decay	SDIRK2	1e-6	53	0.76	0.56
  forced_decay	SDIRK2	1e-8	6340	227	0.60
  forced_decay	Kvaerno3	1e-8	†	†	0.42
  forced_decay	Kvaerno4	1e-8	†	†	0.086
  stiff Prothero–Robinson, $\lambda = -10^4$	TRBDF2	1e-6	0.52	0.00044	0.46
  stiff Prothero–Robinson, $\lambda = -10^4$	TRBDF2	1e-8	2.8	0.0025	0.82
  stiff Prothero–Robinson, $\lambda = -10^4$	SDIRK2	1e-6	19	0.14	184
  stiff Prothero–Robinson, $\lambda = -10^4$	SDIRK2	1e-8	973	0.28	5590
  stiff Prothero–Robinson, $\lambda = -10^4$	Kvaerno3	1e-8	†	†	29
  stiff Prothero–Robinson, $\lambda = -10^4$	Kvaerno4	1e-6	†	†	2.1
  stiff Prothero–Robinson, $\lambda = -10^4$	Kvaerno4	1e-8	†	†	246

  († OrdinaryDiffEq’s Kvaerno3/4 are dominated by their stage-time defect, §5.4, on these time-dependent problems.) On the non-stiff problem $h\gamma J$ is small, the adopted estimate is close to the unsmoothed one, and it removes OrdinaryDiffEq’s loss of accuracy. On the stiff problem it does so for TRBDF2 only: SDIRK2, and Kvaerno3/4 at tight tolerances, still end far above the tolerance, SDIRK2 worse than with OrdinaryDiffEq’s estimate. The unsmoothed estimate meets the tolerance there but takes 5–250 times as many steps (Kvaerno3/4 hit the $10^6$-step limit). Adaptive step sequences differ from OrdinaryDiffEq’s defaults; the oracle checks impose our step sequence, so they are unaffected. Not reported upstream. v1 keeps this estimate for all four methods and recommends TRBDF2 for stiff adaptive leaves; SDIRK2 and Kvaerno3/4 are documented as unreliable there at tight tolerances, and the estimate is revisited at M4b (D59, §13 item 5).

• divided_difference (BackwardEuler; OrdinaryDiffEqSDIRK src/generic_imex_perform_step.jl:2401-2420). With the previous accepted step’s start $(u_{n-1}, t_{n-1})$ in the same leaf invocation:

  e_j = \tfrac{7}{12}\, h^2 \left| \frac{u_{n+1,j} - u_{n,j}}{h\,(t_n + h - t_{n-1})} - \frac{u_{n,j} - u_{n-1,j}}{(t_n - t_{n-1})(t_n + h - t_{n-1})} \right| .
  

  The step’s end time enters as $t_n + h$, as OrdinaryDiffEq writes it, not as the stop time a stop step’s end is set to (§10). In arithmetic: $d_1 = h,((t_n + h) - t_{n-1})$, $d_2 = (t_n - t_{n-1}),((t_n + h) - t_{n-1})$, $r = \tfrac{7}{12},(h \cdot h)$ with $\tfrac{7}{12}$ the table’s divided_difference.constant, and $e_j = r \cdot \big|(u_{n+1,j} - u_{n,j})/d_1 - (u_{n,j} - u_{n-1,j})/d_2\big|$.

  When there is no previous step, the scaled error is set to exactly 1 (accepted). The history resets at every invocation of the leaf, for the reason of §3.3, and after a skipped tiny stop step (§10). After a domain clip (§11) it keeps the accepted state as clipped, which is the next step’s $u_n$. Assumption: resetting per invocation, rather than carrying $u_{n-1}$ across maps and other leaves, is the right reading for split schedules.

6. IMEX Euler and the θ-method

• IMEXEuler: $u_{n+1} = u_n + h f_E(u_n) + h f_I(u_{n+1})$ with $f = f_E + f_I$ by the leaf’s imex partition (run-config §4.5). The implicit equation is solved by §5.2’s Newton with $W = I - hJ_I$, the explicit part evaluated once at $(u_n, t_n)$ and the implicit part at $t_n + h$. As one stage of §5 with $\gamma = 1$: $\tau = u_n + h \cdot f_E(u_n, t_n)$ componentwise, and $\tau = u_n$ when the leaf has no explicit part (no $f_E$ evaluation); $z = h, f_I(\tau + z, t_n + h)$ from $z^{(0)} = 0$; $u_{n+1} = \tau + z$. $W$ uses $J_I$, the Jacobian of the implicit part alone, at $(u_n, t_n)$. With implicit_solve: "linearized" the step is one Rosenbrock–Euler step, $(I - hJ_I(u_n)),\Delta = h,(f_E(u_n) + f_I(u_n))$, $u_{n+1} = u_n + \Delta$ (D43): a one-stage Rosenbrock step with $d = 0$. $f_E$, $f_I$ and $J_I$ are all evaluated at $t_n$, and there is no time-derivative term. In arithmetic, $M = \sigma I - J_I(u_n, t_n)$ with $\sigma = 1/h$ is factored as in §5.3, $\Delta = M^{-1}(f_E + f_I)$ with the two parts evaluated separately and summed componentwise, $f_E$ first, and $u_{n+1} = u_n + \Delta$. Each evaluation of one part, $f_E$ or $f_I$, counts as one right-hand-side call; each part sums its operators in the leaf’s order.
• Theta: $u_{n+1} = u_n + h[(1-\theta) f(u_n, t_n) + \theta f(u_{n+1}, t_n + h)]$ with theta from the leaf (required, run-config §4.2), solved by §5.2’s Newton with $W = I - h\theta J$. θ = 0 is forward Euler (no solve); θ = 1 is backward Euler. As one stage of §5 with $\gamma = \theta$: $\tau = u_n + (h(1-\theta)) \cdot f(u_n, t_n)$, with the coefficient $h(1-\theta)$ rounded once; at θ = 1 the term is skipped and $f(u_n)$ is not evaluated. $z = h, f(\tau + \theta z, t_n + h)$ from $z^{(0)} = 0$, with $J$ at $(u_n, t_n)$, and $u_{n+1} = \tau + \theta z$. At θ = 0, $u_{n+1} = u_n + h \cdot f(u_n, t_n)$.

Neither has an error estimate (D52). See adjoint.md §4.4–§4.5.

7. Error norm

For an error vector $e$ from a step $u_n \to u_{n+1}$ over the leaf’s DOFs $j = 1..N$ (architecture §3.2 pass 3):

\mathrm{sc}_j = \mathrm{abstol}_j + \mathrm{reltol}\cdot\max(|u_{n,j}|, |u_{n+1,j}|), \qquad
E = \sqrt{\frac{1}{N}\sum_{j=1}^{N} \Big(\frac{e_j}{\mathrm{sc}_j}\Big)^2} .

This is DiffEqBase’s calculate_residuals (src/calculate_residuals.jl:9-14) under its default RMS norm (src/common_defaults.jl:69), and KPP’s ros_ErrorNorm.

• Precision (D13, D41). The ratios, squares, sum and square root are computed in f64 whatever the state precision. A Dual number with f64 value type counts as f64; with step-size tangents off its partials are stripped before the norm (adjoint.md §2.3).
• Order (D10). The sum runs over ascending DOF index in the leaf’s plan order, sequentially. Implementations MUST NOT reorder it, use pairwise or SIMD-tree summation, or contract it into fused multiply-adds unless the other implementation does the same; bitwise reproducibility within a language for a fixed decomposition depends on it.
• Scope. With dt_scope: global the sum runs over all the leaf’s DOFs on all ranks, reduced in fixed order (architecture §7.4). With shard or cell it runs over the DOFs of one shard or one cell, and each has its own controller state (§8.6).
• Non-finite. A non-finite $E$ rejects the step, and $h \leftarrow h \cdot q_{\min}$. Deviation: OrdinaryDiffEq has no such branch; a NaN estimate there propagates into $h$ and ends the solve with DtNaN.

8. Step-size control

The leaf’s control.controller picks PI (the default, as in OrdinaryDiffEqCore default_controller, src/alg_utils.jl:740), I or PID. The constants are the controller field of the method’s table: OrdinaryDiffEqCore’s defaults for that algorithm (src/alg_utils.jl:416-889).

Below, $E$ is the step’s scaled error (§7) and $h$ its size. $q$ is the factor that divides $h$ (OrdinaryDiffEq’s convention). $q_{\max}$ means qmax_first_step until the leaf (or cell, under dt_scope: cell) has accepted its first step of the run, and qmax afterwards (src/integrators/controllers.jl:288-293).

8.1 Accept or reject

A step is accepted when $E \le 1$ (PI and I) or when the PID factor is at least accept_safety (§8.4), and also:

• no domain violation (§11), no Newton failure (§5.2), and a non-singular $W$ (§4.2);
• or the step is already at $h \le$ dtmin and the leaf was told to force it (not in v1: there is no force_dtmin option).

8.2 PI controller

OrdinaryDiffEqCore src/integrators/controllers.jl:805-843:

q_{11} = E^{\beta_1}, \qquad q = \mathrm{clamp}\Big(\big(q_{11} / E_{\text{old}}^{\beta_2}\big) / \gamma,\ \frac{1}{q_{\max}},\ \frac{1}{q_{\min}}\Big), \qquad q = 1/q_{\max} \text{ if } E = 0 .

• On acceptance: if $q_{\text{steady,min}} \le q \le q_{\text{steady,max}}$ then $q = 1$; $E_{\text{old}} \leftarrow \max(E, \texttt{qoldinit})$; the proposal is $h/q$.
• On rejection: $h \leftarrow h / \min(1/q_{\min},\ q_{11}/\gamma)$.
• $E_{\text{old}}$ starts at qoldinit.
• Powers are exact (^, powf). Deviation: OrdinaryDiffEqCore uses FastPower.fastpower, which is accurate to about $10^{-4}$.

8.3 I controller

OrdinaryDiffEqCore src/integrators/controllers.jl:666-700: with $k$ = embedded_order,

q = \mathrm{clamp}\big(E^{1/(k+1)}/\gamma,\ 1/q_{\max},\ 1/q_{\min}\big), \qquad q = 1/q_{\max} \text{ if } E = 0 .

On acceptance the dead band applies as in §8.2 and the proposal is $h/q$. On rejection $h \leftarrow h/q$.

8.4 PID controller

Söderlind and Wang (2006), OrdinaryDiffEqCore src/integrators/controllers.jl:1012-1078. With $E’ = \max(E, \epsilon_{f64})$, $\varepsilon_1 = 1/E’$, the previous two accepted steps’ $\varepsilon_2, \varepsilon_3$ (initially 1), and $k = \min(\texttt{order}, \texttt{embedded_order}) + 1$:

\phi = \varepsilon_1^{\beta_1/k}\, \varepsilon_2^{\beta_2/k}\, \varepsilon_3^{\beta_3/k}, \qquad \phi \leftarrow 1 + \arctan(\phi - 1) .

$\beta$ defaults to (beta1, beta2, pid_beta3). Accept when $\phi \ge$ pid_accept_safety. On acceptance: if $q_{\text{steady,min}} \le 1/\phi \le q_{\text{steady,max}}$ then $\phi = 1$; shift $\varepsilon_3 \leftarrow \varepsilon_2$, $\varepsilon_2 \leftarrow \varepsilon_1$; the proposal is $h\phi$. On rejection $h \leftarrow h\phi$. The limiter bounds $\phi$ to about $[0.21, 2.57]$, so qmin and qmax do not apply.

Informative. The default PID betas are the PI betas divided by $k$; that is OrdinaryDiffEq’s default, not a tuned PID.

8.5 Other rejections

• Domain violation (§11) or non-finite $E$ (§7): $h \leftarrow h \cdot q_{\min}$, as OrdinaryDiffEq does for isoutofdomain (OrdinaryDiffEqCore src/integrators/integrator_utils.jl:136, :616).
• Newton failure (§5.2): $h \leftarrow h/\texttt{failfactor}$. The next attempt’s Newton rate starts at 1, as every attempt’s does (§5.2 step 4).
• Singular $W$ (§4.2): always a rejection, with the size the controller’s rejection rule gives at $E = 2$.

The controller’s error history ($E_{\text{old}}$, the PID $\varepsilon$s) changes only on acceptance.

8.6 Proposal, bounds and state

• Proposal. After an accepted step at time $t_{n+1}$ with proposal $\hat h$: $\hat h \leftarrow (t_{n+1} + \hat h) - t_{n+1}$, rounded in f64 (OrdinaryDiffEqCore src/integrators/integrator_utils.jl:1215-1231), then clamped to $[\max(\texttt{dtmin}, \epsilon(t_{n+1})),\ \texttt{dtmax}]$.
• Shortened steps. Let a stop step (§10) be accepted with size $h_s$, shortened from the controller’s size $h$. The controller (§8.2–§8.4) computes the next proposal from $h_s’ = \max(h_s, \texttt{dtmin}, \epsilon(t_n))$, with $t_n$ the step’s start, ($h_s’/q$, or $h_s’\phi$ under PID) when the leaf’s previous step attempt was accepted, and from $h$ when this attempt was the leaf’s first of the run or followed a rejection (D61). This is what OrdinaryDiffEqCore 4.19.0 does: the restore at src/integrators/integrator_utils.jl:633-638 takes the size that modify_dt_for_tstops! saved, and after an accepted step that function runs twice (apply_step!, :191, then loopheader!, :123), so the second call saves the shortened size, after fix_dt_at_bounds! (:122, :1265-1278) has raised it to $\max(\texttt{dtmin}, \epsilon(t))$; on the first step and after a rejection it runs once and the unshortened size is restored. An invocation always ends on an accepted step, so at the start of a later invocation the previous attempt counts as accepted, and no extra state is kept. On van der Pol with Tsit5 at reltol 1e-8 this reproduces OrdinaryDiffEq’s accepted, rejected and right-hand-side counts exactly. Informative. Under splitting every invocation ends on a stop, so a short landing step makes the next invocation start from a short proposal, which then grows by at most qmax per step. That costs extra steps in every macro step (in every cell under dt_scope: cell); D61 accepts the cost for parity with OrdinaryDiffEq.
• State. Each leaf (each shard or cell, under dt_scope: shard or cell) keeps its proposal, $E_{\text{old}}$ or the PID history, and whether it has accepted a step. The Newton rate $\theta_{\text{prev}}$ is not part of it: it lives within one step attempt (§5.2 step 4, D63), and $\eta_{\text{old}}$ has no effect (§5.2 step 5). The state persists across the leaf’s invocations and macro steps, and is part of checkpoints and restarts (architecture §10, §11). A fresh run starts it at the initial values above, with the first step from §9.
• dtmax defaults to the length of the leaf’s current interval; dtmin defaults to 0.

9. First step

When control.dt0 is given, the leaf’s first step of the run is dt0, clamped as in §8.6. Otherwise it is OrdinaryDiffEqCore’s _ode_initdt_oop (src/initdt.jl:354-468), at the leaf’s first invocation, with $p$ = order, and in f64. In it, dtmax is capped at the distance from $t_0$ to the invocation’s first stop (§10), as _determine_initdt passes it (src/integrators/integrator_interface.jl:679-690); dt0 is not capped, since the stop rule shortens it:

1. $\mathrm{sk}_j = \mathrm{abstol}j + |u{0,j}|,\mathrm{reltol}$; $d_0 = \lVert u_0/\mathrm{sk} \rVert$, $f_0 = f(u_0, t_0)$, $d_1 = \lVert f_0/\mathrm{sk} \rVert$ (RMS norms).
2. $\delta_{\min}$ = the next float above $\max(\texttt{dtmin}, \epsilon(t_0))$; $h_{\text{small}} = \max(\delta_{\min}, 10^{-6})$.
3. If $f_0$ has a NaN, or $d_1$ is NaN: return $\delta_{\min}$.
4. $h_0 = h_{\text{small}}$ if $d_0 < 10^{-5}$ or $d_1 < 10^{-5}$, else $h_0 = (d_0/d_1)/100$; then $h_0 \leftarrow \min(h_0, \texttt{dtmax})$.
5. $u_1 = u_0 + h_0 f_0$, $f_1 = f(u_1, t_0 + h_0)$. If $f_1 = f_0$ exactly, return $\max(\delta_{\min}, 100 h_0)$.
6. $d_2 = \lVert (f_1 - f_0)/\mathrm{sk} \rVert / h_0$. If $\max(d_1, d_2) \le 10^{-15}$, $h_1 = \max(h_{\text{small}}, 10^{-3} h_0)$; else $h_1 = 10^{-(2 + \log_{10}\max(d_1, d_2))/p}$.
7. Return $\max(\delta_{\min}, \min(100 h_0, h_1, \texttt{dtmax}))$.

A bind_dt: "leaf" parameter has no step size to take during these evaluations, so an adaptive leaf with one must give control.dt0 (run-config §8 rule 15).

Under dt_scope: cell the first step is computed per cell. Assumption: the two extra $f$ evaluations are per cell too; they are pointwise kernels, so this costs two lane-pool evaluations.

10. Landing on interval ends and tstops

A leaf advances over its node’s interval $[t, t+\Delta]$ (run-config §4). The interval end is a stop, and so is every tstop inside the interval (run-config §5); forcing-epoch boundaries are among them, so no step crosses an epoch (adjoint.md §2.4).

Before each step attempt, with $h$ the controller’s size and $d$ the distance to the next stop (OrdinaryDiffEqCore src/integrators/integrator_utils.jl:271-327):

• $\mathrm{tol} = 100,\epsilon(\max(|t|, |t_{\text{stop}}|))$;
• if $h + \mathrm{tol} < d$, the step is an ordinary one;
• otherwise the step is a stop step of size $\min(h, d)$, and when accepted its end time is set to $t_{\text{stop}}$ exactly, not to $t + h$ (fixed_t_for_tstop_error!, :1032-1038).

Fixed-dt leaves follow the same rule: substeps of dt, the last shortened to land on the stop (run-config §4.2). A stop step shorter than $\epsilon(|t|)$ is accepted without computing it (handle_tstop_step!, :329-337): $t$ becomes $t_{\text{stop}}$ and the step counts as accepted in the statistics, but the proposal, the controller’s error history and whether the leaf has accepted a step are unchanged. The step is recorded on the tape like any accepted step, and replay recognizes it by the same rule. OrdinaryDiffEq instead runs its controller on the previous attempt’s $E$ and the tiny size, which makes the next step tiny; the spec departs from it here (D61).

An ordinary step’s end time is $t_n + h$, accumulated in f64. Over a long fixed-dt run the accumulated time can drift past the landing tolerance, leaving a final step of order $10^{-13}$ (van der Pol at dt 0.01 on $[0, 6.3]$ takes 631 steps, as in OrdinaryDiffEq); both implementations accumulate the same way so step counts agree.

11. Domain enforcement

domain: {nonnegative, on_violation, check_stages} (run-config §4.2, architecture §5.3, D19):

• Check. After a step attempt, and before the accept decision, every listed DOF of $u_{n+1}$ must be $\ge 0$. A negative value is a violation. NaN is not a violation here; it shows up as a non-finite $E$ (§7) or, at fixed dt, as Unstable (§12).
• reject (default): a violation rejects the step (§8.5). A fixed-dt leaf cannot shrink, so a violation there ends the leaf with DomainViolation (§12).
• clip: after acceptance, $u_{n+1,j} \leftarrow \max(u_{n+1,j}, 0)$ on the listed DOFs, and the mask is recorded for the adjoint (adjoint.md §5.7). The clipped state invalidates any stored first-same-as-last evaluation (§3.3).
• check_stages (Rosenbrock only): a negative listed component of any stage state $U_i$ is also a violation (the pre-dawn NO failure, architecture §5.3). A stage violation always rejects, under clip too, since a stage cannot be clipped.

Informative. OrdinaryDiffEq’s isoutofdomain is the same reject-and-shrink rule; clip and check_stages have no OrdinaryDiffEq counterpart.

12. Limits and return codes

A leaf returns one code per invocation (per cell under dt_scope: cell, reduced to the first non-success in cell order). Names follow SciMLBase’s ReturnCode where one fits (SciMLBase src/retcodes.jl; the checks follow check_error, src/integrator_interface.jl:988-1037):

A non-success code from any leaf ends the run with that code, naming the leaf (and cell); the schedule does not retry it.

13. Decisions and assumptions

Decided 2026-10-23 (D55–D59):

1. Kvaerno3/4 stage times (§5.4). The tables use the row sums. Tier-2 comparisons against OrdinaryDiffEq’s Kvaerno methods are limited to autonomous problems. Not reported upstream.
2. KPP’s own step-size rule (§4.4) becomes a controller option, controller: "KPP", at M7 for WRF-Chem parity. The tables already record it.
3. Jacobian reuse in SDIRK (§5.3). v1 forms $J$ on every step attempt. OrdinaryDiffEq-style reuse is revisited at M4b, when eqice’s column solves give real costs.
4. Smoothed SDIRK error estimate (§5.5). Hosea and Shampine’s formula, because OrdinaryDiffEq’s scaled version gives low-accuracy results. Not reported upstream. On stiff problems this fixes TRBDF2 only; see item 5.
5. SDIRK error estimate on stiff problems (§5.5). v1 keeps Hosea and Shampine’s estimate for every SDIRK method. TRBDF2 is the recommended SDIRK method for stiff adaptive leaves. SDIRK2 and Kvaerno3/4 remain available but are unreliable on stiff problems at tight tolerances: on the stiff Prothero–Robinson problem SDIRK2 ended 184–5590 times above the tolerance and Kvaerno3/4 29–246 times at tol = 1e-8, while the unsmoothed estimate met the tolerance at up to 250 times the steps. Revisited at M4b, when eqice’s column solves give real costs; candidates then are smoothing only the stiff components, a different embedded pair for SDIRK2, or an unsmoothed estimate with a step cap. M1’s stiff adaptive cases use TRBDF2, and Kvaerno4 at reltol = 1e-6 (2.1 times the tolerance there).

Assumptions: 6. The SDIRK initial guess (previous stage’s $z$, zero for the first implicit stage) changes converged stages only at the Newton tolerance (§5.2). 7. BackwardEuler’s divided-difference history resets at every invocation of the leaf (§5.5). 8. The first-step algorithm runs per cell under dt_scope: cell (§9).

Decided 2026-10-24 (D63):

9. The Newton rate lives within one step attempt (§5.2 step 4). $\theta_{\text{prev}}$ starts at 1 in every attempt and carries across its stages only, so the step is a function of $(u_n, h)$ and replays bitwise from the tape. OrdinaryDiffEq carries it across steps; the difference is in iteration counts and in converged stages at the Newton tolerance.

Deviations from OrdinaryDiffEq, collected: stage times of Kvaerno3/4 (§5.4); the explicit first stage of ESDIRK methods is evaluated, not reused from $z_s/h$ (§5.1); $u_{n+1} = u_n + \sum_i b_i z_i$ also for stiffly accurate methods (§5.1); the smoothed SDIRK error estimate (§5.5); no Jacobian reuse in SDIRK (§5.3); the Newton rate resets at every step attempt (§5.2, D63); the SDIRK initial guess (§5.2); Rosenbrock23’s error estimate adds $h d, f_t$ (§4.3, D62); a stop step shorter than $\epsilon(|t|)$ leaves the controller unchanged (§10, D61); non-finite error estimates reject (§7); the DomainViolation return code (§12); the grouping of stage sums and embedded estimates, and no fused multiply-adds (§3.1, §3.2, §4.1, §5.1); exact powers in the PI controller (§8.2).