Stiff Kinetics — Read the Solver Report

MODEL "Stiff Kinetics — Read the Solver Report"
DESCRIPTION "Robertson reaction benchmark: compare local tolerance controls, inspect numerical work and check conserved total population."
VERSION "1.0.0"
AUTHOR "Grid Team"
TAGS "chemistry", "kinetics", "ode", "solver-reports"

# Synthetic normalized populations A, B, C; time is in seconds. The constants
# define this benchmark; they are not fitted experimental rate constants.
input EndTime = 100
  VALIDATE BETWEEN 0 AND 100 MESSAGE "Use an endpoint from 0 to 100 seconds for this benchmark."

DEFINE Rates(time, state) = [
  -0.04 * INDEX(state,1,1) + 1e4 * INDEX(state,1,2) * INDEX(state,1,3),
  0.04 * INDEX(state,1,1) - 1e4 * INDEX(state,1,2) * INDEX(state,1,3) - 3e7 * INDEX(state,1,2)^2,
  3e7 * INDEX(state,1,2)^2
]

# These are three deliberate solves. Reading report fields or opening a result
# inspector reuses the stored result; it does not request another integration.
output Automatic = ODE_REPORT((time,state) => Rates(time,state), [1,0,0], 0, EndTime, {format: "record", method: "auto", rtol: 1e-6, atol: 1e-10, max_steps: 100000})
output Stiff = ODE_REPORT((time,state) => Rates(time,state), [1,0,0], 0, EndTime, {format: "record", method: "stiff", rtol: 1e-6, atol: 1e-10, max_steps: 100000})
output Tighter = ODE_REPORT((time,state) => Rates(time,state), [1,0,0], 0, EndTime, {format: "record", method: "stiff", rtol: 1e-10, atol: 1e-14, max_steps: 100000})

output PopulationA = INDEX(Tighter.state,1,1)
output PopulationB = INDEX(Tighter.state,1,2)
output PopulationC = INDEX(Tighter.state,1,3)
output AutomaticConservationResidual = SUM(Automatic.state) - 1
output StiffConservationResidual = SUM(Stiff.state) - 1
output TighterConservationResidual = SUM(Tighter.state) - 1
output LargestToleranceDifference = MAX(ABS(INDEX(Stiff.state,1,1)-PopulationA), ABS(INDEX(Stiff.state,1,2)-PopulationB), ABS(INDEX(Stiff.state,1,3)-PopulationC))
output Interpretation = "A small conservation residual alone does not establish state accuracy. Compare terminal states as tolerances tighten, inspect the actual method and work, and use an independent reference. Tolerance settings control local error estimates; equal settings across methods do not imply equal global errors."
output ReferenceAtHundredSeconds = "At EndTime = 100 only: populations are approximately A = 0.6172348824, B = 0.00000615359127, C = 0.3827589640. Compare values with the retained independent-reference evidence; these rounded numbers are not certified bounds."

END MODEL