EarthSciIntegrator: survey of what the integrator must do

Status: survey, 2026-10-13; §6 (library inventory) added 2026-10-14 from the reuse audit. These are the facts the design rests on; they are not the design itself.

Sources: ../EqWeFiC (WRF, WRF-Chem, CFBM fire), ../eqice (PISM), ../reseact.esm (a CTM that contains the prototype integrator), and the upstream compiler stack (../EarthSciAST, ../EarthSciASTDiff, ../EarthSciASTSplitter.jl, ../EarthSciMLBase, ../EarthSciDiscretizations). Paths are relative to the parent code/ directory.

1. Where each model stands

2. Capability matrix

Legend: ● = needed now, ○ = needed later or optional, blank = not needed.

3. Input contract today (EarthSciAST)

• One flat semi-discrete ODE per document.
  • Julia: f!(du,u,p,t), in place and ForwardDiff-generic.
  • Rust: debug_eval_rhs(state, t, params: &HashMap<String,f64>, force_scalar) plus an allocation-free debug_eval_rhs_into (simulate_array/driver.rs:527-533, :613-621), both #[doc(hidden)]. Its in-crate solve is diffsol 0.11 (Bdf, tr_bdf2, tsit45; driver.rs:2160-2207) with a dense finite-difference Jacobian costing 2n RHS evaluations (:2022-2122), a fresh solver per forcing segment (run_segmented, :1593-1716), and no mass matrix, events, sensitivities or dt_min guard (docs/content/rfcs/rust-diffsol-solver-guard-gap.md).
  • Reactant: an out-of-place (u, p, t, bufs) -> du.
• Layout is variable (species) major: scalars first, then one column-major block per array variable (EarthSciAST/pkg/EarthSciAST.jl/src/tree_walk/state_layout.jl). Per-cell block solves therefore vectorize across cells, struct-of-arrays style; reseact’s blocksolve exploits this.
• Flattening sums operator_compose terms and does not record which operator each term came from.
  • Julia tracks one owner per equation internally (flatten.jl:1050-1056), and Python exposes one source_system per equation (flatten.py:1733-1737). Neither is per-term provenance.
  • Splitting today happens downstream, by locality, in EarthSciASTSplitter.jl (Julia only).
• No Jacobian, sparsity or mass matrix is attached to the compiled problem.
  • EarthSciASTDiff can produce symbolic bands, detect_structure (block_diagonal / banded / general), and per-cell Markowitz LU orderings.
  • None of that is wired into a solver.
• The solver block offers stiffness, abstol/reltol and splitting ∈ {none, lie, strang}, all advisory and document-wide. The RFC explicitly excludes algorithm names, DAE declarations and per-component placement.
• Variable classification includes algebraic_unknowns, but only the :mtk compiler runs them, and :mtk refuses data and PDEs.
• Cadence partition (CONST / DISCRETE / CONTINUOUS) already identifies values that change only at refresh events.
• API conventions to keep:
  • esm_problem → init / step! / solve with SciML spelling;
  • the tolerance chain (call site → document → default);
  • tstops unioned and callback replaced;
  • retcode, name-indexed solutions, remake, ensembles.

4. Other prior art

EarthSciMLBase (MTK-based; does not use .esm)

• SolverStrangThreads / SolverStrangSerial are really first-order Lie splitting with a fixed Δt.
• They run one integrator per cell via reinit!, with static thread chunks and shared mutable parameters (which force a lock in EarthSciData).
• SolverIMEX is a SplitODEProblem with a block-diagonal stiff part.

Worth keeping:

• One compiled per-cell kernel, f(du, u_cell, p, t, c1, c2, c3), with coordinates passed as arguments, run through a pluggable map backend (broadcast, threads, AcceleratedKernels).
• A fused batched LU over (n,n,ncells) blocks (EarthSciMLBase/src/blockdiagonal.jl).
• Data refresh as a preset-time event that also fires at init.

Broken or missing:

• the multi-operator sum (it returns only the last operator’s du);
• OperatorJacobian, vjp, sparse Jacobians, Reactant;
• MPI;
• support for anything but exactly 3 spatial dims.

EarthSciIO

• The provider contract is materialize, refresh(t), refresh_times() and prefetch; refresh_times() becomes the tstops.
• With records_per_sample=2 a provider returns a record bracket, and the model does the time blend.
• prefetch only warms the disk cache, so there is no async overlap.
• Output: Julia has a streaming Zarr writer with a reserved region= argument for subdomain writes. Rust’s Zarr writer writes a whole dataset in one call; a streaming handle is “future work”.

Other

• earthsciml/aq_demo/strang_problem.jl is a hand-rolled Strang split with per-thread Rosenbrock23.
• moves.rs and netcdf-rust contain no time integration.

5. Lessons from the prototype (reseact.esm), as design rules

1. The adjoint is of the realized step sequence. Record the accepted (t, dt) tape and replay with the controller off; re-deriving the sequence “silently takes the adjoint of a different trajectory”.
2. Reverse mode cannot cross a data-dependent trip count. Keep adaptive control on the host, or bound it on device with a static cap and masking.
3. Any parameter that is differentiated must be an explicit argument. A closure over p makes dJ/dp “silently vanish”.
4. Symbolic block Jacobians beat FD and colored AD. They are exact, safe under reverse mode, and 4–25× cheaper to compile. FD has “no good h” across 30 orders of magnitude.
5. Positivity comes from rejection, not modification. The clamp zeroes gradients where it bites.
6. Forcing: refreshes must fire on a set independent of macro_dt, using half-open intervals. Refresh cost dominates at high resolution.
7. A global chemistry dt costs about 6× the per-cell ideal, but every per-cell scheme tried so far lost to per-call overhead at 4×5°. That will probably flip at 0.25°.
8. No caches that grow with grid × run length. Compilation and memory must become per-tile, not per-domain.
9. Treat non-finite values as rejections, and design against XLA:CPU races (a 1-in-2000 NaN) and buffer donation.
10. Upstream SciML pain:
    • similar(jac_prototype) densified BlockDiagonal;
    • the splitting package dropped callbacks and tolerances;
    • a rejected outer step crashed;
    • Reactant cannot trace through the SciML stack.

6. What already exists, verified (reuse audit, 2026-10-14)

Read, not built. Paths are relative to code/; depot paths are under ~/.julia/packages and ~/.cargo/registry/src.

6.1 In the sibling repos

6.2 Julia packages in the depot

6.3 Rust crates in the registry