KPP A Software Environment for the Simulation of Chemical Kinetics - PowerPoint PPT Presentation

KPP – A Software Environment for the Simulation of Chemical Kinetics Adrian Sandu Computer Science Department Virginia Tech NIST, May 25, 2004

Acknowledgements � Dr. Valeriu Damian � Dr. Florian Potra � Dr. Gregory Carmichael � Dr. Dacian Daescu � Dr. Tianfeng Chai � Dr. Mirela Damian Part of this work is supported by � NSF ITR AP&IM-0205198 May 25, 2004

Related Work CHEMKIN ( http://www.reactiondesign.com/lobby/open) Translates symbolic chemical system in a data file that is then used by internal libraries for simulation. Gas-phase kinetics, surface kinetics, reversible equations, transport, mixing, deposition for different types of reactors, direct sensitivity analysis (Senkin). Database of reaction data, graphical postprocessor for results. KINALC ( http://www.chem.leeds.ac.uk/Combustion/kinalc.htm) Postprocessor to CHEMKIN for sensitivity and uncertainty analysis, parameter estimation, and mechanism reduction; etc. KINTECUS ( http://www.kintecus.com) Compiler/ chemical modeling software. Can run heterogeneous and equilibrium chemistry, generates analytical Jacobians, fit/optimize rate constants, initial concentrations, etc. from data; sensitivity analysis; Excel interface. Can use Chemkin models and databases. CANTERA ( http://rayleigh.cds.caltech.edu/~goodwin/cantera) Object-oriented package for chemically-reacting flows. C++ kernel, STL, standard numerical libraries, Interfaces for MATLAB, Python, C++, and Fortran. Capabilities: homogeneous and heterogeneous kinetics, equilibria, reactor modeling, multicomponent transport. LARKIN/LIMKIN ( http://www.zib.de/nowak/codes/limkin_1.0/full) Simulation of LARge systems of chemical reaktion KINetics and parameter identification. Parser generates Fortran code for function and Jacobian, or internal data arrays. DYNAFIT ( http://www.biokin.com/dynafit) Performs nonlinear least-squares regression of chemical kinetic, enzyme kinetic, or ligand- receptor binding data using experimental data. Parses symbolic equations. May 25, 2004

KPP in a Nutshell The Kinetic PreProcessor � Purpose: automatically implement building blocks for large-scale � simulations Parses chemical mechanisms � Generates Fortran and C code for simulation, and direct and adjoint � sensitivity analysis Function, Jacobian, Hessian, Stoichiometric matrix, derivatives of � function&Jacobian w.r.t. rate coefficients Treatment of sparsity � Comprehensive library of numerical integrators � Used in several countries by academia/research/industry � Free! http://www.cs.vt.edu/~asandu/Software/KPP � May 25, 2004

KPP Architecture C C kinetic description files Fortran Fortran Integrator file Utility file atoms species C C equations Fortran Fortran options Driver file I/O file inline code KPP Scanner / Parser Error checking Compute best sparsity structure Reorder species for Substitution best sparsity pre-processor Compute expression trees for assignments C Fortran Code generation C Fortran prod/destr function integrator function utility functions jacobian code driver function(main) I/O functions sparsity structure May 25, 2004

KPP Example SUBROUTINE FunVar ( V, R, F, RCT, A_VAR ) #INCLUDE atoms INCLUDE 'small.h' REAL*8 V(NVAR), R(NRAD), F(NFIX) #DEFVAR REAL*8 RCT(NREACT), A_VAR(NVAR) O = O; O1D = O; O3 = O + O + O; C A - rate for each equation NO = N + O; NO2 = N + O + O; REAL*8 A(NREACT) C Computation of equation rates A(1) = RCT(1)*F(2) #DEFFIX A(2) = RCT(2)*V(2)*F(2) O2 = O + O; M = ignore; A(3) = RCT(3)*V(3) A(4) = RCT(4)*V(2)*V(3) A(5) = RCT(5)*V(3) #EQUATIONS { Small Stratospheric } A(6) = RCT(6)*V(1)*F(1) O2 + hv = 2O : 2.6E-10*SUN**3; A(7) = RCT(7)*V(1)*V(3) O + O2 = O3 : 8.0E-17; A(8) = RCT(8)*V(3)*V(4) O3 + hv = O + O2 : 6.1E-04*SUN; A(9) = RCT(9)*V(2)*V(5) A(10) = RCT(10)*V(5) O + O3 = 2O2 : 1.5E-15; C Aggregate function O3 + hv = O1D + O2 : 1.0E-03*SUN**2; A_VAR(1) = A(5)-A(6)-A(7) O1D + M = O + M : 7.1E-11; A_VAR(2) = 2*A(1)-A(2)+A(3)-A(4)+A(6)- O1D + O3 = 2O2 : 1.2E-10; &A(9)+A(10) A_VAR(3) = A(2)-A(3)-A(4)-A(5)-A(7)-A(8) NO + O3 = NO2 + O2 : 6.0E-15; A_VAR(4) = -A(8)+A(9)+A(10) NO2 + O = NO + O2 : 1.0E-11; A_VAR(5) = A(8)-A(9)-A(10) RETURN NO2 + hv = NO + O : 1.2E-02*SUN; END May 25, 2004

Sparse Jacobians E.g. SAPRC-99 74+5 spc./211 react. IDEAS: NZ=839, NZLU=920 • Chem. interactions: sparsity pattern (off-line) • Min. fill-in reordering 1 • Expand sparsity structure 10 • Row compressed form • Doolitle LU factorization 20 • Loop-free substitution 30 40 #JACOBIAN [ ON | OFF | SPARSE ] 50 JacVar(…), JacVar_SP(…), JacVar_SP_Vec(…), JacVarTR_SP_Vec(…) 60 KppDecomp(…) 70 KppSolve(…), KppSolveTR(…) 74 1 10 20 30 40 50 60 70 74 May 25, 2004

Computational Efficiency (1/Lapack) Dec Sol Dec+7Sol Harwell 0.61 0.21 0.35 Linear Algebra KPP 0.23 0.06 0.12 5 4.5 S−rodas Lsodes rodas 4 Number of Accurate Digits 3.5 3 Stiff I ntegrators 2.5 S−sdirk4 vode sdirk4 2 1.5 S−vode chemeq qssa 1 0.5 8 9 10 20 30 50 80 CPU time [seconds] May 25, 2004

Sparse Hessians E.g. SAPRC99. NZ = 848x2 (0.2%) ∂ 2 f = i H ∂ ∂ i , j , k y y j k • 3-tensor • sparse coordinate format • account for symmetry #HESSIAN [ ON | OFF ] HessVar(…) HessVar_Vec(…) HessVarTR_Vec(…) May 25, 2004

Stoichiometric Form #STOICMAT [ ON | OFF ] STOICM (column compressed) dFunVar_dRcoeff(…) ReactantProd(…) dJacVar_dRcoeff(…) JacVarReactantProd(…) SAPRC−99. Stoichiometric Matrix. NZ = 1024 ( 6 % ) 1 20 Species 40 60 79 1 20 40 60 80 100 120 140 160 180 200 211 Reaction May 25, 2004

Requirements for Numerics � Numerical stability (stiff chemistry) � Accuracy: medium-low (relerr~10 -6 -10 -2 ) � Low Computational Time � Mass Balance � Positivity May 25, 2004

Stiff Integration Methods ( ) − = k k ∑ ∑ α β − [ n ] n i [ n ] n i BDF y h f y i n i = = i 0 i 0 ( ) s ∑ + = + n 1 n y y b f Y , j j Implicit = j 1 Runge-Kutta ( ) s ∑ = + n Y y a f Y i i , j j = j 1 May 25, 2004

Rosenbrock Methods 4 3.5 s ∑ + = + n 1 n y y m k j j = 3 Ros3 j 1 − i 1 ∑ = + n Y y a k 2.5 i i , j j = j 1 SDA ( ) − 2 i 1 ( ) ∑ − ⋅ = + n 1 1 I J k f Y c k Rodas3 γ i i i , j j h h = 1.5 j 1 � No Newton Iterations 1 VODE � Suitable for Stiff Systems TWOSTEP RODAS 0.5 � Mass Conservative SEULEX � Efficient: Low/Med Accuracy 0 0.6 0.8 9 1 2 3 4 5 6 7 8 10 CPU time [seconds] May 25, 2004

Direct Decoupled Sensitivity ( ) ′ = = ∂ ∂ ⎧ y f t , y , p S y p Problem l l ⎨ ( ) ( ) ′ = ⋅ + ≤ ≤ l S J t , y , p S f t , y , p 1 m ⎩ l l p l − γ = ⇒ T I h J P LU ⎡ ⎤ T L P L 0 0 [ ] ⎡ ⎤ ⎢ ⎥ L U 0 ( ) − − γ + 1 T L h JS J U P L 0 ⎢ ⎥ ⎢ ⎥ − γ ℑ = ⋅ 1 p y M O M I h 1 ⎢ ⎥ ⎢ ⎥ M O M 0 [ ] ⎢ ⎥ ⎢ ⎥ L ⎣ ⎦ 0 U ( ) − γ + − 1 T L ⎢ ⎥ h JS J U 0 P L ⎣ ⎦ m p y m May 25, 2004

DDM with KPP − k 1 ∑ + = − + β + n 1 n i n 1 y y h f = i 0 [ ] − k 1 BDF ∑ − β + ⋅ + = − + β + ≤ ≤ n 1 n 1 n i n 1 l I h J S S h f , 1 m l l p (Dunker, ’84) l = i 0 − s s i 1 ∑ ∑ ∑ + = + + = + = + l n 1 n 0 n 1 n i n 0 y y m k , S S m k , Y y a k l l i i i i ij j = = = i 1 i 1 j 1 [ ] ( ) − i 1 ∑ − ⋅ = + + γ n 0 i i 0 n 1 1 I J k f T , Y , p c k h f γ i ij j i t h h = j 1 Rosenbrock ⎛ ⎞ [ ] ( ) ( ) − i 1 ∑ (Sandu et. al. ‘02) ⎜ ⎟ − ⋅ = ⋅ + + l l n i i n i i 1 I J k J T , Y , p S a k f T , Y , p ⎜ ⎟ γ l i ij j p h ⎝ ⎠ l = j 1 − ( ) i 1 ∑ + + ⋅ + × ⋅ + γ + γ l n 0 n n 0 n n n 1 c k J k H S k h J S h f l l ij j p i i i t i p , t h l l = j 1 May 25, 2004

Adjoint Sensitivity Analysis Problem : Stiff ODE, scalar functional. ( ) ( ) ( ) = ≤ ≤ ψ ⇒ ∇ ψ 0 f f y ' f t , y , t t t , y ( t ) y ( t ) 0 y Continuous : Take adjoint of problem, then discretize. ( ) ( ) ( ) ( ) ( ) λ = − ⋅ λ ≥ ≥ λ = ∇ ψ ⇒ ∇ ψ = λ T f 0 f f 0 ' J t , y , t t t , t y ( t ) t y 0 y Discrete : Discretize the problem, then take adjoint: ( ) + = = − k 1 k k L y F y , k 0 , 1 , , N 1 ( ) ( ) ( ) ( ) ( ) λ = ∇ λ = ∇ ψ ⇒ ∇ ψ = λ k k T k N N N 0 F y , y y 0 y y y Note : In both approaches the forward solution needs to be precomputed and stored. May 25, 2004