1 /******************************************************************************
3 KPP - The Kinetic PreProcessor
4 Builds simulation code for chemical kinetic systems
6 Copyright (C) 1995-1996 Valeriu Damian and Adrian Sandu
7 Copyright (C) 1997-2005 Adrian Sandu
9 KPP is free software; you can redistribute it and/or modify it under the
10 terms of the GNU General Public License as published by the Free Software
11 Foundation (http://www.gnu.org/copyleft/gpl.html); either version 2 of the
12 License, or (at your option) any later version.
14 KPP is distributed in the hope that it will be useful, but WITHOUT ANY
15 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
19 You should have received a copy of the GNU General Public License along
20 with this program; if not, consult http://www.gnu.org/copyleft/gpl.html or
21 write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
22 Boston, MA 02111-1307, USA.
25 Computer Science Department
26 Virginia Polytechnic Institute and State University
28 E-mail: sandu@cs.vt.edu
30 ******************************************************************************/
58 ATOM_DEF AtomTable
[ MAX_ATNR
];
59 SPECIES_DEF SpeciesTable
[ MAX_SPECIES
];
60 CODE ReverseCode
[ MAX_SPECIES
];
61 CODE Code
[ MAX_SPECIES
];
67 int Reactive
[ MAX_SPECIES
];
69 INLINE_KEY InlineKeys
[] = { { F77_GLOBAL
, APPEND
, "F77_GLOBAL" },
70 { F77_INIT
, APPEND
, "F77_INIT" },
71 { F77_DATA
, APPEND
, "F77_DATA" },
72 { F77_UTIL
, APPEND
, "F77_UTIL" },
73 { F77_RATES
, APPEND
, "F77_RATES" },
74 { F77_RCONST
, APPEND
, "F77_RCONST" },
75 { F90_GLOBAL
, APPEND
, "F90_GLOBAL" },
76 { F90_INIT
, APPEND
, "F90_INIT" },
77 { F90_DATA
, APPEND
, "F90_DATA" },
78 { F90_UTIL
, APPEND
, "F90_UTIL" },
79 { F90_RATES
, APPEND
, "F90_RATES" },
80 { F90_RCONST
, APPEND
, "F90_RCONST" },
81 { C_GLOBAL
, APPEND
, "C_GLOBAL" },
82 { C_INIT
, APPEND
, "C_INIT" },
83 { C_DATA
, APPEND
, "C_DATA" },
84 { C_UTIL
, APPEND
, "C_UTIL" },
85 { C_RATES
, APPEND
, "C_RATES" },
86 { C_RCONST
, APPEND
, "C_RCONST" },
87 { MATLAB_GLOBAL
, APPEND
, "MATLAB_GLOBAL" },
88 { MATLAB_INIT
, APPEND
, "MATLAB_INIT" },
89 { MATLAB_DATA
, APPEND
, "MATLAB_DATA" },
90 { MATLAB_UTIL
, APPEND
, "MATLAB_UTIL" },
91 { MATLAB_RATES
, APPEND
, "MATLAB_RATES" },
92 { MATLAB_RCONST
, APPEND
, "MATLAB_RCONST" }
96 int useJacobian
= JAC_LU_ROW
;
103 int useDummyindex
= 0;
105 int useLang
= F77_LANG
;
106 int useStochastic
= 0;
108 char integrator
[ MAX_PATH
] = "none";
109 char driver
[ MAX_PATH
] = "none";
110 char runArgs
[ MAX_PATH
] = "";
112 /* mz_rs_20050701+ */
113 /* char varDefault[ MAX_IVAL ] = "1.E-8"; */
114 /* char fixDefault[ MAX_IVAL ] = "1.E-8"; */
115 /* double cfactor = 1.09E+10; */
116 char varDefault
[ MAX_IVAL
] = "0.";
117 char fixDefault
[ MAX_IVAL
] = "0.";
119 /* mz_rs_20050701- */
121 ATOM crtAtoms
[ MAX_ATOMS
];
124 char *fileList
[ MAX_FILES
];
127 double Abs( double x
)
129 return x
> 0 ? x
: -x
;
132 void DefineInitializeNbr( char *cmd
)
136 n
= sscanf( cmd
, "%d", &initNr
);
138 ScanError("Bad number of species to initialize <%s>", cmd
);
141 void DefineXGrid( char *cmd
)
146 n
= sscanf( cmd
, "%d", &xNr
);
148 ScanError("Bad X grid number <%s>", cmd
);
151 void DefineYGrid( char *cmd
)
156 n
= sscanf( cmd
, "%d", &yNr
);
158 ScanError("Bad Y grid number <%s>", cmd
);
161 void DefineZGrid( char *cmd
)
166 n
= sscanf( cmd
, "%d", &zNr
);
168 ScanError("Bad Z grid number <%s>", cmd
);
171 void CmdFunction( char *cmd
)
173 if( EqNoCase( cmd
, "AGGREGATE" ) ) {
177 if( EqNoCase( cmd
, "SPLIT" ) ) {
181 ScanError("'%s': Unknown parameter for #FUNCTION [AGGREGATE|SPLIT]", cmd
);
184 void CmdJacobian( char *cmd
)
186 if( EqNoCase( cmd
, "OFF" ) ) {
187 useJacobian
= JAC_OFF
;
191 if( EqNoCase( cmd
, "FULL" ) ) {
192 useJacobian
= JAC_FULL
;
196 if( EqNoCase( cmd
, "SPARSE_LU_ROW" ) ) {
197 useJacobian
= JAC_LU_ROW
;
201 if( EqNoCase( cmd
, "SPARSE_ROW" ) ) {
202 useJacobian
= JAC_ROW
;
206 ScanError("'%s': Unknown parameter for #JACOBIAN [OFF|FULL|SPARSE_LU_ROW|SPARSE_ROW]", cmd
);
209 void SparseData( char *cmd
) {
210 ScanError("Deprecated use of #SPARSEDATA %s: see #JACOBIAN for equivalent functionality", cmd
);
213 void CmdHessian( char *cmd
)
215 if( EqNoCase( cmd
, "OFF" ) ) {
219 if( EqNoCase( cmd
, "ON" ) ) {
223 ScanError("'%s': Unknown parameter for #HESSIAN [ON|OFF]", cmd
);
226 void CmdStoicmat( char *cmd
)
228 if( EqNoCase( cmd
, "OFF" ) ) {
232 if( EqNoCase( cmd
, "ON" ) ) {
236 ScanError("'%s': Unknown parameter for #STOICMAT [ON|OFF]", cmd
);
239 void CmdDouble( char *cmd
)
241 if( EqNoCase( cmd
, "OFF" ) ) {
245 if( EqNoCase( cmd
, "ON" ) ) {
249 ScanError("'%s': Unknown parameter for #DOUBLE [ON|OFF]", cmd
);
252 void CmdReorder( char *cmd
)
254 if( EqNoCase( cmd
, "OFF" ) ) {
258 if( EqNoCase( cmd
, "ON" ) ) {
262 ScanError("'%s': Unknown parameter for #REORDER [ON|OFF]", cmd
);
265 void CmdMex( char *cmd
)
267 if( EqNoCase( cmd
, "OFF" ) ) {
271 if( EqNoCase( cmd
, "ON" ) ) {
275 ScanError("'%s': Unknown parameter for #MEX [ON|OFF]", cmd
);
278 void CmdDummyindex( char *cmd
)
280 if( EqNoCase( cmd
, "OFF" ) ) {
284 if( EqNoCase( cmd
, "ON" ) ) {
288 ScanError("'%s': Unknown parameter for #DUMMYINDEX [ON|OFF]", cmd
);
291 void CmdEqntags( char *cmd
)
293 if( EqNoCase( cmd
, "OFF" ) ) {
297 if( EqNoCase( cmd
, "ON" ) ) {
301 ScanError("'%s': Unknown parameter for #EQNTAGS [ON|OFF]", cmd
);
304 void CmdUse( char *cmd
)
306 ScanError("Deprecated command '#USE %s';\nReplace with '#LANGUAGE %s'.",cmd
,cmd
);
310 void CmdLanguage( char *cmd
)
312 if( EqNoCase( cmd
, "FORTRAN77" ) ) {
316 if( EqNoCase( cmd
, "FORTRAN" ) ) {
317 ScanWarning("Fortran version not specified in '#LANGUAGE %s'. Will use Fortran 77.", cmd
);
321 if( EqNoCase( cmd
, "FORTRAN90" ) ) {
325 if( EqNoCase( cmd
, "MATLAB" ) ) {
326 useLang
= MATLAB_LANG
;
329 if( EqNoCase( cmd
, "C" ) ) {
333 ScanError("'%s': Unknown parameter for #LANGUAGE [Fortran77|Fortran90|C|Matlab]", cmd
);
336 void CmdStochastic( char *cmd
)
338 if( EqNoCase( cmd
, "ON" ) ) {
342 if( EqNoCase( cmd
, "OFF" ) ) {
346 ScanError("'%s': Unknown parameter for #STOCHASTIC [OFF|ON]", cmd
);
349 void CmdIntegrator( char *cmd
)
351 strcpy( integrator
, cmd
);
354 void CmdDriver( char *cmd
)
356 strcpy( driver
, cmd
);
359 void CmdRun( char *cmd
)
361 strcpy( runArgs
, cmd
);
364 int FindAtom( char *atname
)
368 for( i
=0; i
<AtomNr
; i
++ )
369 if( EqNoCase( AtomTable
[ i
].name
, atname
) ) {
375 void DeclareAtom( char *atname
)
379 code
= FindAtom( atname
);
381 ScanError("Multiple declaration for atom %s.", atname
);
384 if( AtomNr
>= MAX_ATNR
) {
385 Error("Too many atoms");
389 strcpy( AtomTable
[ AtomNr
].name
, atname
);
390 AtomTable
[ AtomNr
].check
= NO_CHECK
;
391 AtomTable
[ AtomNr
].masscheck
= 0;
395 void SetAtomType( char *atname
, int type
)
399 code
= FindAtom( atname
);
401 ScanError("Undefined atom %s.", atname
);
404 AtomTable
[ code
].check
= type
;
411 for( i
=0; i
<AtomNr
; i
++ ) {
412 if( AtomTable
[ i
].check
!= CANCEL_CHECK
)
413 AtomTable
[ i
].check
= DO_CHECK
;
415 SetAtomType( "IGNORE", NO_CHECK
);
418 void AddAtom( char *atname
, char *nr
)
422 code
= FindAtom( atname
);
424 ScanError("Undefined atom %s.", atname
);
427 crtAtoms
[ crtAtomNr
].code
= (unsigned char)code
;
428 crtAtoms
[ crtAtomNr
].nr
= (unsigned char)atoi(nr
);
432 int FindSpecies( char *spname
)
436 for( i
=0; i
<SpeciesNr
; i
++ )
437 if( EqNoCase( SpeciesTable
[ i
].name
, spname
) ) {
441 if( EqNoCase( SpeciesTable
[ MAX_SPECIES
-1 - i
].name
, spname
) ) {
442 return MAX_SPECIES
-1 - i
;
447 void StoreSpecies( int index
, int type
, char *spname
)
451 strcpy( SpeciesTable
[ index
].name
, spname
);
452 SpeciesTable
[ index
].type
= type
;
453 *SpeciesTable
[ index
].ival
= '\0';
454 SpeciesTable
[ index
].lookat
= 0;
455 SpeciesTable
[ index
].moni
= 0;
456 SpeciesTable
[ index
].trans
= 0;
457 if( (SpeciesTable
[ index
].nratoms
== 0) || ( crtAtomNr
> 0 ) ) {
458 SpeciesTable
[ index
].nratoms
= crtAtomNr
;
459 for( i
= 0; i
< crtAtomNr
; i
++ )
460 SpeciesTable
[ index
].atoms
[i
] = crtAtoms
[i
];
465 void DeclareSpecies( int type
, char *spname
)
469 code
= FindSpecies( spname
);
471 ScanError("Multiple declaration for species %s.", spname
);
474 if( SpeciesNr
>= MAX_SPECIES
) {
475 Error("Too many species");
478 StoreSpecies( SpeciesNr
, type
, spname
);
482 void SetSpcType( int type
, char *spname
)
487 if( EqNoCase( spname
, "VAR_SPEC" ) ) {
488 for( i
= 0; i
< SpeciesNr
; i
++ )
489 if( SpeciesTable
[i
].type
== VAR_SPC
)
490 SpeciesTable
[i
].type
= type
;
493 if( EqNoCase( spname
, "FIX_SPEC" ) ) {
494 for( i
= 0; i
< SpeciesNr
; i
++ )
495 if( SpeciesTable
[i
].type
== FIX_SPC
)
496 SpeciesTable
[i
].type
= type
;
499 if( EqNoCase( spname
, "ALL_SPEC" ) ) {
500 for( i
= 0; i
< SpeciesNr
; i
++ )
501 SpeciesTable
[i
].type
= type
;
505 code
= FindSpecies( spname
);
507 ScanError("Undefined species %s.", spname
);
510 SpeciesTable
[ code
].type
= type
;
513 void AssignInitialValue( char *spname
, char *spval
)
518 if( EqNoCase( spname
, "CFACTOR" ) ) {
519 code
= sscanf( spval
, "%lg", &cf
);
521 ScanWarning("Invalid CFACTOR value: %s", spval
);
528 if( EqNoCase( spname
, "VAR_SPEC" ) ) {
529 strcpy( varDefault
, spval
);
534 if( EqNoCase( spname
, "FIX_SPEC" ) ) {
535 strcpy( fixDefault
, spval
);
539 if( EqNoCase( spname
, "ALL_SPEC" ) ) {
540 strcpy( varDefault
, spval
);
541 strcpy( fixDefault
, spval
);
545 code
= FindSpecies( spname
);
547 ScanError("Undefined species %s.", spname
);
550 strcpy( SpeciesTable
[ code
].ival
, spval
);
553 void StoreEquationRate( char *rate
, char *label
)
560 kreact
= &kr
[ EqnNr
];
561 strcpy( kreact
->label
, label
);
563 kreact
->type
= PHOTO
;
564 strcpy( kreact
->val
.st
, rate
);
568 n
= sscanf( rate
, "%lf%s", &f
, buf
);
570 kreact
->type
= NUMBER
;
574 kreact
->type
= EXPRESION
;
575 strcpy( kreact
->val
.st
, rate
);
584 float atcnt
[ MAX_ATNR
];
590 if( EqnNr
>= MAX_EQN
) {
591 Error("Too many equations");
595 for( i
= 0; i
< AtomNr
; i
++ )
598 for( spc
= 0; spc
< SpcNr
; spc
++ ) {
599 sp
= &SpeciesTable
[ Code
[spc
] ];
600 if( Stoich_Left
[spc
][EqnNr
] != 0 ) {
601 for( i
= 0; i
< sp
->nratoms
; i
++ )
602 atcnt
[ sp
->atoms
[i
].code
] += Stoich_Left
[spc
][EqnNr
] * sp
->atoms
[i
].nr
;
604 if( Stoich_Right
[spc
][EqnNr
] != 0 ) {
605 for( i
= 0; i
< sp
->nratoms
; i
++ )
606 atcnt
[ sp
->atoms
[i
].code
] -= Stoich_Right
[spc
][EqnNr
] * sp
->atoms
[i
].nr
;
613 for( i
= 0; i
< AtomNr
; i
++ ) {
614 if ( Abs( atcnt
[i
] ) > 1e-5 ) {
615 if ( AtomTable
[i
].check
== CANCEL_CHECK
) {
619 if ( AtomTable
[i
].check
== NO_CHECK
) {
622 if ( AtomTable
[i
].check
== DO_CHECK
) {
624 sprintf(errmsg
, "%s %s", errmsg
, AtomTable
[i
].name
);
631 ScanWarning( "(eqn %d) Atom balance mismatch for:%s.", EqnNr
+1, errmsg
);
633 for( j
= 0; j
< SpcNr
; j
++ )
634 if( Stoich_Left
[j
][EqnNr
] != 0 )
635 { index
= j
; break; }
636 for( i
= 0; i
< EqnNr
; i
++ ) {
638 r1
= Stoich_Left
[index
][EqnNr
];
639 r2
= Stoich_Left
[index
][i
];
640 for( j
= 0; j
< SpcNr
; j
++ ) {
641 if( r1
* Stoich_Left
[j
][i
] != r2
* Stoich_Left
[j
][EqnNr
] )
642 { equal
= 0; break; }
643 if( r1
* Stoich_Right
[j
][i
] != r2
* Stoich_Right
[j
][EqnNr
] )
644 { equal
= 0; break; }
648 ScanError( "Duplicate equation: "
649 " (eqn<%d> = eqn<%d> )", i
+1, EqnNr
+1 );
651 ScanError( "Linearly dependent equations: "
652 "( %.0f eqn<%d> = %.0f eqn<%d> )",
653 r1
, i
+1, r2
, EqnNr
+1 );
660 void ProcessTerm( int side
, char *sign
, char *coef
, char *spname
)
668 code
= FindSpecies( spname
);
670 ScanError("Undefined species %s.", spname
);
674 crtSpec
= ReverseCode
[ code
];
676 if(EqNoCase(spname
,"HV")) isPhoto
= 1;
678 if ( crtSpec
== NO_CODE
) {
679 if( MAX_SPECIES
- code
<= 2 ) falseSpcNr
++;
681 Code
[ crtSpec
] = code
;
682 ReverseCode
[ code
] = crtSpec
;
687 sscanf( buf
, "%lf", &val
);
690 case LHS
: Stoich_Left
[ crtSpec
][ EqnNr
] += val
;
691 Stoich
[ crtSpec
][ EqnNr
] -= val
;
692 Reactive
[ crtSpec
] = 1;
694 case RHS
: Stoich_Right
[ crtSpec
][ EqnNr
] += val
;
695 Stoich
[ crtSpec
][ EqnNr
] += val
;
700 void AddLumpSpecies( char *spname
)
704 code
= FindSpecies( spname
);
706 ScanError("Undefined species %s.", spname
);
714 void CheckLump( char *spname
)
718 code
= FindSpecies( spname
);
720 ScanError("Undefined species %s.", spname
);
728 void AddLookAt( char *spname
)
732 code
= FindSpecies( spname
);
734 ScanError("Undefined species %s.", spname
);
738 SpeciesTable
[ code
].lookat
= 1;
745 for( i
=0; i
<SpeciesNr
; i
++ )
746 SpeciesTable
[ i
].lookat
= 1;
749 void AddMonitor( char *spname
)
753 code
= FindSpecies( spname
);
755 SpeciesTable
[ code
].moni
= 1;
759 code
= FindAtom( spname
);
761 AtomTable
[ code
].masscheck
= 1;
765 ScanError("Undefined species or atom %s.", spname
);
768 void AddTransport( char *spname
)
772 code
= FindSpecies( spname
);
774 ScanError("Undefined species %s.", spname
);
778 SpeciesTable
[ code
].trans
= 1;
785 for( i
=0; i
<SpeciesNr
; i
++ )
786 SpeciesTable
[ i
].trans
= 1;
789 void AddUseFile( char *fname
)
791 fileList
[fileNr
] = (char*)malloc(strlen(fname
)+1);
792 strcpy(fileList
[fileNr
], fname
);
796 char * AppendString( char * s1
, char * s2
, int * maxlen
, int addlen
)
803 s1
= (char*)malloc( *maxlen
);
807 if( strlen( s1
) + strlen( s2
) >= *maxlen
) {
808 s1
= (char*)realloc( (void*)s1
, *maxlen
);
814 char * ReplaceString( char * s1
, char * s2
, int * maxlen
, int addlen
)
820 *maxlen
= strlen( s2
);
821 s1
= (char*)malloc( 1+*maxlen
);
827 void AddInlineCode( char * ctx
, char * s
)
831 int totallength
; /* mz_rs_20050607 */
835 for( i
= 0; i
< INLINE_OPT
; i
++ )
836 if( EqNoCase( ctx
, InlineKeys
[i
].kname
) ) {
837 key
= InlineKeys
[i
].key
;
838 c
= &InlineCode
[key
];
839 type
= InlineKeys
[i
].type
;
843 printf( "\n'%s': Unknown inline option (ignored)", ctx
);
847 /* mz_rs_20050607+ */
849 totallength
= strlen( c
->code
)+strlen( s
);
851 totallength
= strlen( s
);
852 if (totallength
>MAX_INLINE
)
853 ScanError("\nInline code for %s is too long (%d>%d).\nIncrease MAX_INLINE in scan.h and recompile kpp!",
854 ctx
, totallength
, MAX_INLINE
);
855 /* mz_rs_20050607- */
858 case APPEND
: c
->code
= AppendString( c
->code
, s
, &c
->maxlen
, MAX_INLINE
);
860 case REPLACE
: c
->code
= ReplaceString( c
->code
, s
, &c
->maxlen
, MAX_INLINE
);
865 int ParseEquationFile( char * filename
)
870 for( i
= 0; i
< MAX_SPECIES
; i
++ ) {
871 ReverseCode
[i
] = NO_CODE
;
874 for( i
= 0; i
< MAX_SPECIES
; i
++ ) {
875 for( j
= 0; j
< MAX_EQN
; j
++ ) {
876 Stoich_Left
[i
][j
] = 0;
878 Stoich_Right
[i
][j
] = 0;
881 for( i
= 0; i
< MAX_SPECIES
; i
++ ) {
882 SpeciesTable
[ i
].nratoms
= 0;
885 for( i
= 0; i
< INLINE_OPT
; i
++ ) {
886 InlineCode
[i
].code
= NULL
;
887 InlineCode
[i
].maxlen
= 0;
893 DeclareAtom( "CANCEL" );
894 SetAtomType( "CANCEL", CANCEL_CHECK
);
895 DeclareAtom( "IGNORE" );
896 SetAtomType( "IGNORE", NO_CHECK
);
897 DeclareSpecies( DUMMY_SPC
, "???" );
898 StoreSpecies( MAX_SPECIES
-1, DUMMY_SPC
, "HV" );
899 AddAtom( "CANCEL", "1" );
900 StoreSpecies( MAX_SPECIES
-2, DUMMY_SPC
, "PROD" );
902 code
= Parser( filename
);