VODE numerical Jacobian corrupts later columns when abundances are renormalized

Suggested priority: Medium
Reviewed revision: 3ac0bda4f00da3eed0446beadfc1cd5c51b99a40

Problem and impact

The finite-difference loop perturbs one component and restores only that component after rhs(). With integrator.renormalize_abundances = true, clean_state()` changes every species component during each RHS call.

Consequently, subsequent columns are differenced from a modified composition while still subtracting the original savf. The stored integration state also changes during Jacobian evaluation. This is an order-one derivative error, not just a small composition-roundoff effect. It applies to numerical Jacobian selection, including a retry that switches to numerical differentiation.

Affected code

• integration/VODE/vode_dvjac.H:120
• integration/integrator_type_strang.H:70

Reproduction and observed behavior

From the repository root, run:

python3 agent-notes/integration-review-2026-09-06/reproduce.py

Run vode_jac. Start at (X1, X2, e) = (0.5, 0.5, 10) and use an RHS that returns the state after the real Strang clean_state() has normalized the species. The exact derivative of its first component is (0.5, -0.5, 0).

y after=0.5000000009313226,0.4999999990686774 expected=.5,.5
J(1,2)=-0.2500000018626451 expected=-.5; J(1,3)=0.00625 expected=0

The test leaves LU factorization as a no-op so that the matrix assembled by dvjac() can be inspected directly; it does not test the linear solver.

The companion script and captured output are part of this review. These are focused CPU tests with minimal AMReX/physics substitutes, not full network or GPU integration tests. Attach the script when submitting this issue.

Suggested fix and regression coverage

Save the entire baseline solution after the baseline RHS call. Restore all components before and after every perturbation, as the Rosenbrock numerical-Jacobian implementation already does. Verify every column against the derivative of the normalized RHS, and verify that the baseline state is preserved.

Start with integration/VODE/vode_dvjac.H:120 and compare its numerical-Jacobian state handling with the Rosenbrock implementation; also inspect clean_state() at integration/integrator_type_strang.H:70. Run python3 agent-notes/integration-review-2026-09-06/reproduce.py and inspect the assembled matrix. Done means every column matches the normalized RHS derivative and the baseline state remains unchanged.