[AI prototype] PositiveIntegrators.jl

[visr]

This is a little test balloon to check if https://github.com/NumericalMathematics/PositiveIntegrators.jl is worth considering.
It is on top of both #3039 and #3082, relevant commit is ad72990. Still untested, and PID is not yet supported. AI analysis:

What PositiveIntegrators.jl offers

Modified Patankar-Runge-Kutta (MPRK) schemes that unconditionally guarantee positivity of solutions by requiring the ODE to be expressed as a Production-Destruction System (PDS):

$$u_i'(t) = \sum_j \left( p_{ij}(t,u) - d_{ij}(t,u) \right), \quad p_{ij} \geq 0,\ d_{ij} \geq 0$$

Available solvers: MPE (1st order), MPRK22 (2nd, adaptive), MPRK43I/MPRK43II (3rd, adaptive), SSPMPRK22/SSPMPRK43, MPDeC (arbitrary order). All are linearly implicit and integrate with the SciML ecosystem.

Current Ribasim approach to positivity

1. isoutofdomain (util.jl): rejects any step where u.basin[i] < 0, forcing dt reduction
2. Clamping (solve.jl): ifelse(u.basin[i] > 0, u.basin[i], 0.0) when computing levels/areas
3. Low-storage reduction factor: smoothly reduces evaporation/infiltration as storage → 0
4. Stiff implicit solvers (QNDF, Rodas5P) with auto-sparse Jacobians via DifferentiationInterface.jl

Natural fit with Ribasim's architecture

Ribasim's apply_flow_to_basins! in solve.jl already has PDS semantics:

function apply_flow_to_basins!(du, q, inflow_id, outflow_id)
    if inflow_id.type == NodeType.Basin
        du.basin[inflow_id.idx] -= q   # destruction from source
    end
    if outflow_id.type == NodeType.Basin
        du.basin[outflow_id.idx] += q   # production into destination
    end
end

This is precisely a production term $p_{\text{dest},\text{src}} = q$ when $q > 0$. The water balance is a sum of such pairwise flows plus external sources/sinks — the exact structure a PDS requires.

Sparse matrix support — not a bottleneck

The p_prototype keyword in PDSProblem accepts SparseMatrixCSC. From the linear advection benchmark, sparse + KLUFactorization() achieved a 260x speedup over dense (45.8s → 174ms for N=1000, near-zero allocations).

For Ribasim this maps cleanly:

• Sparsity pattern = network adjacency (known at construction from the graph topology)
• Nonzeros = $O(\text{links})$, not $O(N_{\text{basins}}^2)$
• KLU is already a dependency — Ribasim uses KLUFactorization for its implicit solvers
• Bidirectional links (e.g. LinearResistance) need entries in both directions in the sparsity pattern

Reformulation effort

Production matrix P!(P, u, p, t) — Moderate

Each flow computation currently ends with apply_flow_to_basins!(du, q, src, dst). The PDS equivalent:

Current code	PDS equivalent
du.basin[dst.idx] += q (q > 0)	P[dst.idx, src.idx] = q
du.basin[src.idx] -= q (q > 0)	(implied by off-diagonal convention)
du.basin[dst.idx] += q (q < 0, bidirectional)	P[src.idx, dst.idx] = -q
Precipitation into basin i	P[i, i] += precip (diagonal source)
Flow boundary into basin i	P[i, i] += inflow (diagonal source)

Destruction vector D!(d, u, p, t) — Simple

External sinks that leave the system entirely:

• d[i] += evaporation[i]
• d[i] += infiltration[i]

Sign-splitting for bidirectional flows — Straightforward

LinearResistance and ManningResistance can flow in either direction based on head difference. At each evaluation:

if q >= 0
    P[outflow_basin_idx, inflow_basin_idx] = q
else
    P[inflow_basin_idx, outflow_basin_idx] = -q
end

The sparsity pattern must include both directions for these links.

Key challenges

Challenge	Severity	Mitigation
PID integral states	Medium	These aren't volumes and can be negative. Options: (a) use std_rhs fallback for non-MPRK solvers on the full system while MPRK only sees basin states; (b) operator splitting; (c) shift integral states by an offset to keep them positive (fragile). Best option is likely (a): restructure so PID integrals are handled outside the PDS, e.g. as a separate coupled equation.
Bidirectional flow sign-splitting	Low	Runtime conditional in P! fill; sparsity pattern includes both directions. Cheap.
std_rhs for fallback	Low	Already have water_balance! — pass it as the std_rhs keyword so non-MPRK algorithms (QNDF, Rodas5P) still work unchanged.
Callbacks & discrete control	Low	PDSProblem <: AbstractODEProblem — all SciML callbacks work.
Adaptive time stepping	None	MPRK22, MPRK43I/II, MPDeC all support adaptive dt.
Custom state vector type	Medium	Ribasim uses CArray (named-tuple-indexed vector). Need to verify PositiveIntegrators.jl doesn't assume plain Vector. The p_prototype is only for the P matrix; u should work if it's <: AbstractVector.

What you gain

1. Guaranteed positivity — no more step rejections from isoutofdomain, no more ifelse clamping
2. Remove workarounds — the low_storage_threshold reduction factor becomes unnecessary for positivity (may still be physically motivated)
3. Better mass conservation — MPRK methods are designed to preserve $\sum u_i$ for conservative systems; even for non-conservative PDS, the production-destruction bookkeeping reduces numerical mass balance errors
4. Simpler solver diagnostics — no "negative storage" errors to debug

What you lose / trade off

1. Automatic sparse Jacobian via DifferentiationInterface.jl — replaced by manual P! matrix fill (but this is arguably more transparent)
2. High-order stiff solvers (QNDF, Rodas5P) — MPRK methods top out at 3rd order (or arbitrary via MPDeC). For very stiff systems the implicit BDF methods may still be faster per accuracy. However, std_rhs lets you keep them as options.
3. Some refactoring — splitting the water_balance! logic into production-matrix fill vs. standard RHS

Implementation path

Phase 1: Proof of concept (basin-only, no PID)

1. Add PositiveIntegrators to Project.toml
2. Build p_prototype sparse matrix from network graph at model construction
3. Implement ribasim_P!(P, u, p, t) and ribasim_D!(d, u, p, t) by refactoring the flow formulation logic
4. Create PDSProblem(...; p_prototype, std_rhs = water_balance!)
5. Test on basic model (no PID) with MPRK22(1.0; linsolve = KLUFactorization())
6. Validate: mass balance, no negative storages, same results as current approach

Phase 2: Handle PID integral states

• Option A: Extend state vector so PID integrals are "phantom" components with $p_{ii} = d_{ii} = 0$ and handle their dynamics via std_rhs (only MPRK for basin components). Requires checking if PositiveIntegrators.jl supports partial PDS (likely not out of the box).
• Option B: Operator splitting — advance basins with MPRK, advance PID integrals with a standard method per step.
• Option C: Add PID integral states to the PDS with a large positive offset ($u_{\text{integral}} = u_{\text{raw}} + C$) so they stay positive. Fragile but simple.

Phase 3: Performance tuning

• Benchmark against current QNDF/Rodas5P on medium-sized models
• Try MPRK43I / MPDeC for higher accuracy
• Compare step counts (MPRK never rejects for positivity, but may need small steps for accuracy)

Phase 4: Full integration

• Add MPRK algorithms to solver config options
• Remove or gate isoutofdomain when MPRK is selected
• Update documentation

Verdict

The fit is natural. Ribasim's water balance is already structured as pairwise flows between basins + external sources/sinks — exactly a PDS. Sparse matrix support with KLU eliminates the performance concern for large models. The main engineering work is:

1. Refactoring water_balance! into a P!/D! matrix-fill form (~moderate effort, same logic restructured)
2. Deciding how to handle PID integral states (~design decision)
3. Validating numerical accuracy matches or improves on current solvers

[gdalle]

Hey there @visr, can you explain why the automatic sparse Jacobian capabilities are lost here? Anything I can do to help?

[visr]

Thanks for the offer @gdalle! I didn't dive into it too much yet, but I think the PDS formulation and algorithms just don't use a Jacobian, but use the production matrix P and D instead to form the linear system.

It is just a little test. For future reference, in the first version the basic testmodel is about 1000x slower. The FunctionCallingAffect seems to suddenly dominate, here is a profile: PDS-basic.html