1 | /****************************************************************************** |
---|
2 | |
---|
3 | KPP - The Kinetic PreProcessor |
---|
4 | Builds simulation code for chemical kinetic systems |
---|
5 | |
---|
6 | Copyright (C) -2021 996 Valeriu Damian and Adrian Sandu |
---|
7 | Copyright (C) -2021 005 Adrian Sandu |
---|
8 | |
---|
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. |
---|
13 | |
---|
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 |
---|
17 | details. |
---|
18 | |
---|
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. |
---|
23 | |
---|
24 | Adrian Sandu |
---|
25 | Computer Science Department |
---|
26 | Virginia Polytechnic Institute and State University |
---|
27 | Blacksburg, VA 24060 |
---|
28 | E-mail: sandu@cs.vt.edu |
---|
29 | |
---|
30 | ******************************************************************************/ |
---|
31 | |
---|
32 | |
---|
33 | #include "gdata.h" |
---|
34 | #include "scan.h" |
---|
35 | #include "y.tab.h" |
---|
36 | #include <stdlib.h> |
---|
37 | #include <string.h> |
---|
38 | #include <math.h> |
---|
39 | |
---|
40 | int AtomNr = 0; |
---|
41 | int SpeciesNr = 0; |
---|
42 | int EqnNr = 0; |
---|
43 | int SpcNr = 0; |
---|
44 | int VarNr = 0; |
---|
45 | int VarActiveNr = 0; |
---|
46 | int FixNr = 0; |
---|
47 | int VarStartNr = 0; |
---|
48 | int FixStartNr = 0; |
---|
49 | |
---|
50 | |
---|
51 | int initNr = -1; |
---|
52 | int xNr = 0; |
---|
53 | int yNr = 0; |
---|
54 | int zNr = 0; |
---|
55 | |
---|
56 | int falseSpcNr = 0; |
---|
57 | |
---|
58 | ATOM_DEF AtomTable[ MAX_ATNR ]; |
---|
59 | SPECIES_DEF SpeciesTable[ MAX_SPECIES ]; |
---|
60 | CODE ReverseCode[ MAX_SPECIES ]; |
---|
61 | CODE Code[ MAX_SPECIES ]; |
---|
62 | KREACT kr[ MAX_EQN ]; |
---|
63 | |
---|
64 | float** Stoich_Left; |
---|
65 | float** Stoich; |
---|
66 | float** Stoich_Right; |
---|
67 | int Reactive[ MAX_SPECIES ]; |
---|
68 | |
---|
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" } |
---|
93 | }; |
---|
94 | |
---|
95 | int useAggregate = 1; |
---|
96 | int useJacobian = JAC_LU_ROW; |
---|
97 | int useJacSparse = 1; |
---|
98 | int useHessian = 1; |
---|
99 | int useStoicmat = 1; |
---|
100 | int useDouble = 1; |
---|
101 | int useReorder = 1; |
---|
102 | int useMex = 1; |
---|
103 | int useDummyindex = 0; |
---|
104 | int useEqntags = 0; |
---|
105 | int useLang = F77_LANG; |
---|
106 | int useStochastic = 0; |
---|
107 | /* if useValues=1 KPP replaces parameters like NVAR etc. |
---|
108 | by their values in vector/matrix declarations */ |
---|
109 | int useDeclareValues = 0; |
---|
110 | |
---|
111 | char integrator[ MAX_PATH ] = "none"; |
---|
112 | char driver[ MAX_PATH ] = "none"; |
---|
113 | char runArgs[ MAX_PATH ] = ""; |
---|
114 | |
---|
115 | /* mz_rs_20050701+ */ |
---|
116 | /* char varDefault[ MAX_IVAL ] = "1.E-8"; */ |
---|
117 | /* char fixDefault[ MAX_IVAL ] = "1.E-8"; */ |
---|
118 | /* double cfactor = 1.09E+10; */ |
---|
119 | char varDefault[ MAX_IVAL ] = "0."; |
---|
120 | char fixDefault[ MAX_IVAL ] = "0."; |
---|
121 | double cfactor = 1.; |
---|
122 | /* mz_rs_20050701- */ |
---|
123 | |
---|
124 | ATOM crtAtoms[ MAX_ATOMS ]; |
---|
125 | int crtAtomNr = 0; |
---|
126 | |
---|
127 | char *fileList[ MAX_FILES ]; |
---|
128 | int fileNr = 0; |
---|
129 | |
---|
130 | double Abs( double x ) |
---|
131 | { |
---|
132 | return x > 0 ? x : -x; |
---|
133 | } |
---|
134 | |
---|
135 | void DefineInitializeNbr( char *cmd ) |
---|
136 | { |
---|
137 | int n; |
---|
138 | |
---|
139 | n = sscanf( cmd, "%d", &initNr); |
---|
140 | if( n != 1 ) |
---|
141 | ScanError("Bad number of species to initialize <%s>", cmd); |
---|
142 | } |
---|
143 | |
---|
144 | void DefineXGrid( char *cmd ) |
---|
145 | { |
---|
146 | int n; |
---|
147 | |
---|
148 | xNr = 1; |
---|
149 | n = sscanf( cmd, "%d", &xNr); |
---|
150 | if( n != 1 ) |
---|
151 | ScanError("Bad X grid number <%s>", cmd); |
---|
152 | } |
---|
153 | |
---|
154 | void DefineYGrid( char *cmd ) |
---|
155 | { |
---|
156 | int n; |
---|
157 | |
---|
158 | yNr = 1; |
---|
159 | n = sscanf( cmd, "%d", &yNr); |
---|
160 | if( n != 1 ) |
---|
161 | ScanError("Bad Y grid number <%s>", cmd); |
---|
162 | } |
---|
163 | |
---|
164 | void DefineZGrid( char *cmd ) |
---|
165 | { |
---|
166 | int n; |
---|
167 | |
---|
168 | zNr = 1; |
---|
169 | n = sscanf( cmd, "%d", &zNr); |
---|
170 | if( n != 1 ) |
---|
171 | ScanError("Bad Z grid number <%s>", cmd); |
---|
172 | } |
---|
173 | |
---|
174 | void CmdFunction( char *cmd ) |
---|
175 | { |
---|
176 | if( EqNoCase( cmd, "AGGREGATE" ) ) { |
---|
177 | useAggregate = 1; |
---|
178 | return; |
---|
179 | } |
---|
180 | if( EqNoCase( cmd, "SPLIT" ) ) { |
---|
181 | useAggregate = 0; |
---|
182 | return; |
---|
183 | } |
---|
184 | ScanError("'%s': Unknown parameter for #FUNCTION [AGGREGATE|SPLIT]", cmd ); |
---|
185 | } |
---|
186 | |
---|
187 | void CmdDeclareValues( char *cmd ) |
---|
188 | { |
---|
189 | if( EqNoCase( cmd, "SYMBOL" ) ) { |
---|
190 | useDeclareValues = 0; |
---|
191 | return; |
---|
192 | } |
---|
193 | if( EqNoCase( cmd, "VALUE" ) ) { |
---|
194 | useDeclareValues = 1; |
---|
195 | return; |
---|
196 | } |
---|
197 | ScanError("'%s': Unknown parameter for #DECLARE [SYMBOL|VALUE]", cmd ); |
---|
198 | } |
---|
199 | |
---|
200 | void CmdJacobian( char *cmd ) |
---|
201 | { |
---|
202 | if( EqNoCase( cmd, "OFF" ) ) { |
---|
203 | useJacobian = JAC_OFF; |
---|
204 | useJacSparse = 0; |
---|
205 | return; |
---|
206 | } |
---|
207 | if( EqNoCase( cmd, "FULL" ) ) { |
---|
208 | useJacobian = JAC_FULL; |
---|
209 | useJacSparse = 0; |
---|
210 | return; |
---|
211 | } |
---|
212 | if( EqNoCase( cmd, "SPARSE_LU_ROW" ) ) { |
---|
213 | useJacobian = JAC_LU_ROW; |
---|
214 | useJacSparse = 1; |
---|
215 | return; |
---|
216 | } |
---|
217 | if( EqNoCase( cmd, "SPARSE_ROW" ) ) { |
---|
218 | useJacobian = JAC_ROW; |
---|
219 | useJacSparse = 1; |
---|
220 | return; |
---|
221 | } |
---|
222 | ScanError("'%s': Unknown parameter for #JACOBIAN [OFF|FULL|SPARSE_LU_ROW|SPARSE_ROW]", cmd ); |
---|
223 | } |
---|
224 | |
---|
225 | void SparseData( char *cmd ) { |
---|
226 | ScanError("Deprecated use of #SPARSEDATA %s: see #JACOBIAN for equivalent functionality", cmd ); |
---|
227 | } |
---|
228 | |
---|
229 | void CmdHessian( char *cmd ) |
---|
230 | { |
---|
231 | if( EqNoCase( cmd, "OFF" ) ) { |
---|
232 | useHessian = 0; |
---|
233 | return; |
---|
234 | } |
---|
235 | if( EqNoCase( cmd, "ON" ) ) { |
---|
236 | useHessian = 1; |
---|
237 | return; |
---|
238 | } |
---|
239 | ScanError("'%s': Unknown parameter for #HESSIAN [ON|OFF]", cmd ); |
---|
240 | } |
---|
241 | |
---|
242 | void CmdStoicmat( char *cmd ) |
---|
243 | { |
---|
244 | if( EqNoCase( cmd, "OFF" ) ) { |
---|
245 | useStoicmat = 0; |
---|
246 | return; |
---|
247 | } |
---|
248 | if( EqNoCase( cmd, "ON" ) ) { |
---|
249 | useStoicmat = 1; |
---|
250 | return; |
---|
251 | } |
---|
252 | ScanError("'%s': Unknown parameter for #STOICMAT [ON|OFF]", cmd ); |
---|
253 | } |
---|
254 | |
---|
255 | void CmdDouble( char *cmd ) |
---|
256 | { |
---|
257 | if( EqNoCase( cmd, "OFF" ) ) { |
---|
258 | useDouble = 0; |
---|
259 | return; |
---|
260 | } |
---|
261 | if( EqNoCase( cmd, "ON" ) ) { |
---|
262 | useDouble = 1; |
---|
263 | return; |
---|
264 | } |
---|
265 | ScanError("'%s': Unknown parameter for #DOUBLE [ON|OFF]", cmd ); |
---|
266 | } |
---|
267 | |
---|
268 | void CmdReorder( char *cmd ) |
---|
269 | { |
---|
270 | if( EqNoCase( cmd, "OFF" ) ) { |
---|
271 | useReorder = 0; |
---|
272 | return; |
---|
273 | } |
---|
274 | if( EqNoCase( cmd, "ON" ) ) { |
---|
275 | useReorder = 1; |
---|
276 | return; |
---|
277 | } |
---|
278 | ScanError("'%s': Unknown parameter for #REORDER [ON|OFF]", cmd ); |
---|
279 | } |
---|
280 | |
---|
281 | void CmdMex( char *cmd ) |
---|
282 | { |
---|
283 | if( EqNoCase( cmd, "OFF" ) ) { |
---|
284 | useMex = 0; |
---|
285 | return; |
---|
286 | } |
---|
287 | if( EqNoCase( cmd, "ON" ) ) { |
---|
288 | useMex = 1; |
---|
289 | return; |
---|
290 | } |
---|
291 | ScanError("'%s': Unknown parameter for #MEX [ON|OFF]", cmd ); |
---|
292 | } |
---|
293 | |
---|
294 | void CmdDummyindex( char *cmd ) |
---|
295 | { |
---|
296 | if( EqNoCase( cmd, "OFF" ) ) { |
---|
297 | useDummyindex = 0; |
---|
298 | return; |
---|
299 | } |
---|
300 | if( EqNoCase( cmd, "ON" ) ) { |
---|
301 | useDummyindex = 1; |
---|
302 | return; |
---|
303 | } |
---|
304 | ScanError("'%s': Unknown parameter for #DUMMYINDEX [ON|OFF]", cmd ); |
---|
305 | } |
---|
306 | |
---|
307 | void CmdEqntags( char *cmd ) |
---|
308 | { |
---|
309 | if( EqNoCase( cmd, "OFF" ) ) { |
---|
310 | useEqntags = 0; |
---|
311 | return; |
---|
312 | } |
---|
313 | if( EqNoCase( cmd, "ON" ) ) { |
---|
314 | useEqntags = 1; |
---|
315 | return; |
---|
316 | } |
---|
317 | ScanError("'%s': Unknown parameter for #EQNTAGS [ON|OFF]", cmd ); |
---|
318 | } |
---|
319 | |
---|
320 | void CmdUse( char *cmd ) |
---|
321 | { |
---|
322 | ScanError("Deprecated command '#USE %s';\nReplace with '#LANGUAGE %s'.",cmd,cmd ); |
---|
323 | } |
---|
324 | |
---|
325 | |
---|
326 | void CmdLanguage( char *cmd ) |
---|
327 | { |
---|
328 | if( EqNoCase( cmd, "FORTRAN77" ) ) { |
---|
329 | useLang = F77_LANG; |
---|
330 | return; |
---|
331 | } |
---|
332 | if( EqNoCase( cmd, "FORTRAN" ) ) { |
---|
333 | ScanWarning("Fortran version not specified in '#LANGUAGE %s'. Will use Fortran 77.", cmd); |
---|
334 | useLang = F77_LANG; |
---|
335 | return; |
---|
336 | } |
---|
337 | if( EqNoCase( cmd, "FORTRAN90" ) ) { |
---|
338 | useLang = F90_LANG; |
---|
339 | return; |
---|
340 | } |
---|
341 | if( EqNoCase( cmd, "MATLAB" ) ) { |
---|
342 | useLang = MATLAB_LANG; |
---|
343 | return; |
---|
344 | } |
---|
345 | if( EqNoCase( cmd, "C" ) ) { |
---|
346 | useLang = C_LANG; |
---|
347 | return; |
---|
348 | } |
---|
349 | ScanError("'%s': Unknown parameter for #LANGUAGE [Fortran77|Fortran90|C|Matlab]", cmd ); |
---|
350 | } |
---|
351 | |
---|
352 | void CmdStochastic( char *cmd ) |
---|
353 | { |
---|
354 | if( EqNoCase( cmd, "ON" ) ) { |
---|
355 | useStochastic = 1; |
---|
356 | return; |
---|
357 | } |
---|
358 | if( EqNoCase( cmd, "OFF" ) ) { |
---|
359 | useStochastic = 0; |
---|
360 | return; |
---|
361 | } |
---|
362 | ScanError("'%s': Unknown parameter for #STOCHASTIC [OFF|ON]", cmd ); |
---|
363 | } |
---|
364 | |
---|
365 | void CmdIntegrator( char *cmd ) |
---|
366 | { |
---|
367 | strcpy( integrator, cmd ); |
---|
368 | } |
---|
369 | |
---|
370 | void CmdDriver( char *cmd ) |
---|
371 | { |
---|
372 | strcpy( driver, cmd ); |
---|
373 | } |
---|
374 | |
---|
375 | void CmdRun( char *cmd ) |
---|
376 | { |
---|
377 | strcpy( runArgs, cmd ); |
---|
378 | } |
---|
379 | |
---|
380 | int FindAtom( char *atname ) |
---|
381 | { |
---|
382 | int i; |
---|
383 | |
---|
384 | for( i=0; i<AtomNr; i++ ) |
---|
385 | if( EqNoCase( AtomTable[ i ].name, atname ) ) { |
---|
386 | return i; |
---|
387 | } |
---|
388 | return -1; |
---|
389 | } |
---|
390 | |
---|
391 | void DeclareAtom( char *atname ) |
---|
392 | { |
---|
393 | int code; |
---|
394 | |
---|
395 | code = FindAtom( atname ); |
---|
396 | if ( code >= 0 ) { |
---|
397 | ScanError("Multiple declaration for atom %s.", atname ); |
---|
398 | return; |
---|
399 | } |
---|
400 | if( AtomNr >= MAX_ATNR ) { |
---|
401 | Error("Too many atoms"); |
---|
402 | return; |
---|
403 | } |
---|
404 | |
---|
405 | strcpy( AtomTable[ AtomNr ].name, atname ); |
---|
406 | AtomTable[ AtomNr ].check = NO_CHECK; |
---|
407 | AtomTable[ AtomNr ].masscheck = 0; |
---|
408 | AtomNr++; |
---|
409 | } |
---|
410 | |
---|
411 | void SetAtomType( char *atname, int type ) |
---|
412 | { |
---|
413 | int code; |
---|
414 | |
---|
415 | code = FindAtom( atname ); |
---|
416 | if ( code < 0 ) { |
---|
417 | ScanError("Undefined atom %s.", atname ); |
---|
418 | return; |
---|
419 | } |
---|
420 | AtomTable[ code ].check = type; |
---|
421 | } |
---|
422 | |
---|
423 | void CheckAll() |
---|
424 | { |
---|
425 | int i; |
---|
426 | |
---|
427 | for( i=0; i<AtomNr; i++ ) { |
---|
428 | if( AtomTable[ i ].check != CANCEL_CHECK ) |
---|
429 | AtomTable[ i ].check = DO_CHECK; |
---|
430 | } |
---|
431 | SetAtomType( "IGNORE", NO_CHECK ); |
---|
432 | } |
---|
433 | |
---|
434 | void AddAtom( char *atname, char *nr ) |
---|
435 | { |
---|
436 | int code; |
---|
437 | |
---|
438 | code = FindAtom( atname ); |
---|
439 | if ( code < 0 ) { |
---|
440 | ScanError("Undefined atom %s.", atname ); |
---|
441 | return; |
---|
442 | } |
---|
443 | crtAtoms[ crtAtomNr ].code = (unsigned char)code; |
---|
444 | crtAtoms[ crtAtomNr ].nr = (unsigned char)atoi(nr); |
---|
445 | crtAtomNr++; |
---|
446 | } |
---|
447 | |
---|
448 | int FindSpecies( char *spname ) |
---|
449 | { |
---|
450 | int i; |
---|
451 | |
---|
452 | for( i=0; i<SpeciesNr; i++ ) |
---|
453 | if( EqNoCase( SpeciesTable[ i ].name, spname ) ) { |
---|
454 | return i; |
---|
455 | } |
---|
456 | for( i=0; i<2; i++ ) |
---|
457 | if( EqNoCase( SpeciesTable[ MAX_SPECIES -1 - i ].name, spname ) ) { |
---|
458 | return MAX_SPECIES -1 - i; |
---|
459 | } |
---|
460 | return -1; |
---|
461 | } |
---|
462 | |
---|
463 | void StoreSpecies( int index, int type, char *spname ) |
---|
464 | { |
---|
465 | int i; |
---|
466 | |
---|
467 | strcpy( SpeciesTable[ index ].name, spname ); |
---|
468 | SpeciesTable[ index ].type = type; |
---|
469 | *SpeciesTable[ index ].ival = '\0'; |
---|
470 | SpeciesTable[ index ].lookat = 0; |
---|
471 | SpeciesTable[ index ].moni = 0; |
---|
472 | SpeciesTable[ index ].trans = 0; |
---|
473 | if( (SpeciesTable[ index ].nratoms == 0) || ( crtAtomNr > 0 ) ) { |
---|
474 | SpeciesTable[ index ].nratoms = crtAtomNr; |
---|
475 | for( i = 0; i < crtAtomNr; i++ ) |
---|
476 | SpeciesTable[ index ].atoms[i] = crtAtoms[i]; |
---|
477 | } |
---|
478 | crtAtomNr = 0; |
---|
479 | } |
---|
480 | |
---|
481 | void DeclareSpecies( int type, char *spname ) |
---|
482 | { |
---|
483 | int code; |
---|
484 | |
---|
485 | code = FindSpecies( spname ); |
---|
486 | if ( code >= 0 ) { |
---|
487 | ScanError("Multiple declaration for species %s.", spname ); |
---|
488 | return; |
---|
489 | } |
---|
490 | if( SpeciesNr >= MAX_SPECIES ) { |
---|
491 | Error("Too many species"); |
---|
492 | return; |
---|
493 | } |
---|
494 | StoreSpecies( SpeciesNr, type, spname ); |
---|
495 | SpeciesNr++; |
---|
496 | } |
---|
497 | |
---|
498 | void SetSpcType( int type, char *spname ) |
---|
499 | { |
---|
500 | int code; |
---|
501 | int i; |
---|
502 | |
---|
503 | if( EqNoCase( spname, "VAR_SPEC" ) ) { |
---|
504 | for( i = 0; i < SpeciesNr; i++ ) |
---|
505 | if( SpeciesTable[i].type == VAR_SPC ) |
---|
506 | SpeciesTable[i].type = type; |
---|
507 | return; |
---|
508 | } |
---|
509 | if( EqNoCase( spname, "FIX_SPEC" ) ) { |
---|
510 | for( i = 0; i < SpeciesNr; i++ ) |
---|
511 | if( SpeciesTable[i].type == FIX_SPC ) |
---|
512 | SpeciesTable[i].type = type; |
---|
513 | return; |
---|
514 | } |
---|
515 | if( EqNoCase( spname, "ALL_SPEC" ) ) { |
---|
516 | for( i = 0; i < SpeciesNr; i++ ) |
---|
517 | SpeciesTable[i].type = type; |
---|
518 | return; |
---|
519 | } |
---|
520 | |
---|
521 | code = FindSpecies( spname ); |
---|
522 | if ( code < 0 ) { |
---|
523 | ScanError("Undefined species %s.", spname ); |
---|
524 | return; |
---|
525 | } |
---|
526 | SpeciesTable[ code ].type = type; |
---|
527 | } |
---|
528 | |
---|
529 | void AssignInitialValue( char *spname , char *spval ) |
---|
530 | { |
---|
531 | int code; |
---|
532 | double cf; |
---|
533 | |
---|
534 | if( EqNoCase( spname, "CFACTOR" ) ) { |
---|
535 | code = sscanf( spval, "%lg", &cf ); |
---|
536 | if( code != 1 ) { |
---|
537 | ScanWarning("Invalid CFACTOR value: %s", spval); |
---|
538 | return; |
---|
539 | } |
---|
540 | cfactor = cf; |
---|
541 | return; |
---|
542 | } |
---|
543 | |
---|
544 | if( EqNoCase( spname, "VAR_SPEC" ) ) { |
---|
545 | strcpy( varDefault, spval ); |
---|
546 | return; |
---|
547 | } |
---|
548 | |
---|
549 | |
---|
550 | if( EqNoCase( spname, "FIX_SPEC" ) ) { |
---|
551 | strcpy( fixDefault, spval ); |
---|
552 | return; |
---|
553 | } |
---|
554 | |
---|
555 | if( EqNoCase( spname, "ALL_SPEC" ) ) { |
---|
556 | strcpy( varDefault, spval ); |
---|
557 | strcpy( fixDefault, spval ); |
---|
558 | return; |
---|
559 | } |
---|
560 | |
---|
561 | code = FindSpecies( spname ); |
---|
562 | if ( code < 0 ) { |
---|
563 | ScanError("Undefined species %s.", spname ); |
---|
564 | return; |
---|
565 | } |
---|
566 | strcpy( SpeciesTable[ code ].ival, spval ); |
---|
567 | } |
---|
568 | |
---|
569 | void StoreEquationRate( char *rate, char *label ) |
---|
570 | { |
---|
571 | double f; |
---|
572 | char buf[ MAX_K ]; |
---|
573 | int n; |
---|
574 | KREACT *kreact; |
---|
575 | |
---|
576 | kreact = &kr[ EqnNr ]; |
---|
577 | strcpy( kreact->label, label ); |
---|
578 | if( isPhoto ) { |
---|
579 | kreact->type = PHOTO; |
---|
580 | strcpy( kreact->val.st, rate ); |
---|
581 | isPhoto = 0; |
---|
582 | return; |
---|
583 | } |
---|
584 | n = sscanf( rate, "%lf%s", &f, buf ); |
---|
585 | if ( n == 1 ) { |
---|
586 | kreact->type = NUMBER; |
---|
587 | kreact->val.f = f; |
---|
588 | return; |
---|
589 | } |
---|
590 | kreact->type = EXPRESION; |
---|
591 | strcpy( kreact->val.st, rate ); |
---|
592 | return; |
---|
593 | } |
---|
594 | |
---|
595 | void CheckEquation() |
---|
596 | { |
---|
597 | int i,j; |
---|
598 | int equal, index; |
---|
599 | double r1, r2; |
---|
600 | float atcnt[ MAX_ATNR ]; |
---|
601 | int spc; |
---|
602 | SPECIES_DEF *sp; |
---|
603 | char errmsg[80]; |
---|
604 | int err; |
---|
605 | |
---|
606 | if( EqnNr >= MAX_EQN ) { |
---|
607 | Error("Too many equations"); |
---|
608 | return; |
---|
609 | } |
---|
610 | |
---|
611 | for( i = 0; i < AtomNr; i++ ) |
---|
612 | atcnt[i] = 0; |
---|
613 | |
---|
614 | for( spc = 0; spc < SpcNr; spc++ ) { |
---|
615 | sp = &SpeciesTable[ Code[spc] ]; |
---|
616 | if( Stoich_Left[spc][EqnNr] != 0 ) { |
---|
617 | for( i = 0; i < sp->nratoms; i++ ) |
---|
618 | atcnt[ sp->atoms[i].code ] += Stoich_Left[spc][EqnNr] * sp->atoms[i].nr; |
---|
619 | } |
---|
620 | if( Stoich_Right[spc][EqnNr] != 0 ) { |
---|
621 | for( i = 0; i < sp->nratoms; i++ ) |
---|
622 | atcnt[ sp->atoms[i].code ] -= Stoich_Right[spc][EqnNr] * sp->atoms[i].nr; |
---|
623 | } |
---|
624 | } |
---|
625 | |
---|
626 | *errmsg = 0; |
---|
627 | err = 0; |
---|
628 | |
---|
629 | for( i = 0; i < AtomNr; i++ ) { |
---|
630 | if ( Abs( atcnt[i] ) > 1e-5 ) { |
---|
631 | if ( AtomTable[i].check == CANCEL_CHECK ) { |
---|
632 | err = 0; |
---|
633 | break; |
---|
634 | } |
---|
635 | if ( AtomTable[i].check == NO_CHECK ) { |
---|
636 | continue; |
---|
637 | } |
---|
638 | if ( AtomTable[i].check == DO_CHECK ) { |
---|
639 | err = 1; |
---|
640 | sprintf(errmsg, "%s %s", errmsg, AtomTable[i].name ); |
---|
641 | continue; |
---|
642 | } |
---|
643 | } |
---|
644 | } |
---|
645 | |
---|
646 | if ( err ) |
---|
647 | ScanWarning( "(eqn %d) Atom balance mismatch for:%s.", EqnNr+1, errmsg ); |
---|
648 | |
---|
649 | for( j = 0; j < SpcNr; j++ ) |
---|
650 | if( Stoich_Left[j][EqnNr] != 0 ) |
---|
651 | { index = j; break; } |
---|
652 | for( i = 0; i < EqnNr; i++ ) { |
---|
653 | equal = 1; |
---|
654 | r1 = Stoich_Left[index][EqnNr]; |
---|
655 | r2 = Stoich_Left[index][i]; |
---|
656 | for( j = 0; j < SpcNr; j++ ) { |
---|
657 | if( r1 * Stoich_Left[j][i] != r2 * Stoich_Left[j][EqnNr] ) |
---|
658 | { equal = 0; break; } |
---|
659 | if( r1 * Stoich_Right[j][i] != r2 * Stoich_Right[j][EqnNr] ) |
---|
660 | { equal = 0; break; } |
---|
661 | } |
---|
662 | if ( equal ) { |
---|
663 | if( r1 == r2 ) |
---|
664 | ScanError( "Duplicate equation: " |
---|
665 | " (eqn<%d> = eqn<%d> )", i+1, EqnNr+1 ); |
---|
666 | else |
---|
667 | ScanError( "Linearly dependent equations: " |
---|
668 | "( %.0f eqn<%d> = %.0f eqn<%d> )", |
---|
669 | r1, i+1, r2, EqnNr+1 ); |
---|
670 | break; |
---|
671 | } |
---|
672 | } |
---|
673 | EqnNr++; |
---|
674 | } |
---|
675 | |
---|
676 | void ProcessTerm( int side, char *sign, char *coef, char *spname ) |
---|
677 | { |
---|
678 | int code; |
---|
679 | CODE crtSpec; |
---|
680 | double val; |
---|
681 | char buf[40]; |
---|
682 | |
---|
683 | |
---|
684 | code = FindSpecies( spname ); |
---|
685 | if ( code < 0 ) { |
---|
686 | ScanError("Undefined species %s.", spname ); |
---|
687 | return; |
---|
688 | } |
---|
689 | |
---|
690 | crtSpec = ReverseCode[ code ]; |
---|
691 | |
---|
692 | if(EqNoCase(spname,"HV")) isPhoto = 1; |
---|
693 | |
---|
694 | if ( crtSpec == NO_CODE ) { |
---|
695 | if( MAX_SPECIES - code <= 2 ) falseSpcNr++; |
---|
696 | crtSpec = SpcNr++; |
---|
697 | Code[ crtSpec ] = code; |
---|
698 | ReverseCode[ code ] = crtSpec; |
---|
699 | } |
---|
700 | |
---|
701 | strcpy( buf, sign ); |
---|
702 | strcat( buf, coef ); |
---|
703 | sscanf( buf, "%lf", &val ); |
---|
704 | |
---|
705 | switch( side ) { |
---|
706 | case LHS: Stoich_Left[ crtSpec ][ EqnNr ] += val; |
---|
707 | Stoich[ crtSpec ][ EqnNr ] -= val; |
---|
708 | Reactive[ crtSpec ] = 1; |
---|
709 | break; |
---|
710 | case RHS: Stoich_Right[ crtSpec ][ EqnNr ] += val; |
---|
711 | Stoich[ crtSpec ][ EqnNr ] += val; |
---|
712 | break; |
---|
713 | } |
---|
714 | } |
---|
715 | |
---|
716 | void AddLumpSpecies( char *spname ) |
---|
717 | { |
---|
718 | int code; |
---|
719 | |
---|
720 | code = FindSpecies( spname ); |
---|
721 | if ( code < 0 ) { |
---|
722 | ScanError("Undefined species %s.", spname ); |
---|
723 | return; |
---|
724 | } |
---|
725 | |
---|
726 | /* ... */ |
---|
727 | |
---|
728 | } |
---|
729 | |
---|
730 | void CheckLump( char *spname ) |
---|
731 | { |
---|
732 | int code; |
---|
733 | |
---|
734 | code = FindSpecies( spname ); |
---|
735 | if ( code < 0 ) { |
---|
736 | ScanError("Undefined species %s.", spname ); |
---|
737 | return; |
---|
738 | } |
---|
739 | |
---|
740 | /* ... */ |
---|
741 | |
---|
742 | } |
---|
743 | |
---|
744 | void AddLookAt( char *spname ) |
---|
745 | { |
---|
746 | int code; |
---|
747 | |
---|
748 | code = FindSpecies( spname ); |
---|
749 | if ( code < 0 ) { |
---|
750 | ScanError("Undefined species %s.", spname ); |
---|
751 | return; |
---|
752 | } |
---|
753 | |
---|
754 | SpeciesTable[ code ].lookat = 1; |
---|
755 | } |
---|
756 | |
---|
757 | void LookAtAll() |
---|
758 | { |
---|
759 | int i; |
---|
760 | |
---|
761 | for( i=0; i<SpeciesNr; i++ ) |
---|
762 | SpeciesTable[ i ].lookat = 1; |
---|
763 | } |
---|
764 | |
---|
765 | void AddMonitor( char *spname ) |
---|
766 | { |
---|
767 | int code; |
---|
768 | |
---|
769 | code = FindSpecies( spname ); |
---|
770 | if ( code >= 0 ) { |
---|
771 | SpeciesTable[ code ].moni = 1; |
---|
772 | return; |
---|
773 | } |
---|
774 | |
---|
775 | code = FindAtom( spname ); |
---|
776 | if ( code >= 0 ) { |
---|
777 | AtomTable[ code ].masscheck = 1; |
---|
778 | return; |
---|
779 | } |
---|
780 | |
---|
781 | ScanError("Undefined species or atom %s.", spname ); |
---|
782 | } |
---|
783 | |
---|
784 | void AddTransport( char *spname ) |
---|
785 | { |
---|
786 | int code; |
---|
787 | |
---|
788 | code = FindSpecies( spname ); |
---|
789 | if ( code < 0 ) { |
---|
790 | ScanError("Undefined species %s.", spname ); |
---|
791 | return; |
---|
792 | } |
---|
793 | |
---|
794 | SpeciesTable[ code ].trans = 1; |
---|
795 | } |
---|
796 | |
---|
797 | void TransportAll() |
---|
798 | { |
---|
799 | int i; |
---|
800 | |
---|
801 | for( i=0; i<SpeciesNr; i++ ) |
---|
802 | SpeciesTable[ i ].trans = 1; |
---|
803 | } |
---|
804 | |
---|
805 | void AddUseFile( char *fname ) |
---|
806 | { |
---|
807 | fileList[fileNr] = (char*)malloc(strlen(fname)+1); |
---|
808 | strcpy(fileList[fileNr], fname); |
---|
809 | fileNr++; |
---|
810 | } |
---|
811 | |
---|
812 | char * AppendString( char * s1, char * s2, int * maxlen, int addlen ) |
---|
813 | { |
---|
814 | char * tmp; |
---|
815 | |
---|
816 | *maxlen += addlen; |
---|
817 | |
---|
818 | if( !s1 ) { |
---|
819 | s1 = (char*)malloc( *maxlen ); |
---|
820 | *s1 = 0; |
---|
821 | } |
---|
822 | |
---|
823 | if( strlen( s1 ) + strlen( s2 ) >= *maxlen ) { |
---|
824 | s1 = (char*)realloc( (void*)s1, *maxlen ); |
---|
825 | } |
---|
826 | strcat( s1, s2 ); |
---|
827 | return s1; |
---|
828 | } |
---|
829 | |
---|
830 | char * ReplaceString( char * s1, char * s2, int * maxlen, int addlen ) |
---|
831 | { |
---|
832 | char * tmp; |
---|
833 | |
---|
834 | if( s1 ) free(s1); |
---|
835 | |
---|
836 | *maxlen = strlen( s2 ); |
---|
837 | s1 = (char*)malloc( 1+*maxlen ); |
---|
838 | strcpy( s1, s2 ); |
---|
839 | |
---|
840 | return s1; |
---|
841 | } |
---|
842 | |
---|
843 | void AddInlineCode( char * ctx, char * s ) |
---|
844 | { |
---|
845 | ICODE * c; |
---|
846 | int i, key, type; |
---|
847 | int totallength; /* mz_rs_20050607 */ |
---|
848 | |
---|
849 | c = NULL; |
---|
850 | |
---|
851 | for( i = 0; i < INLINE_OPT; i++ ) |
---|
852 | if( EqNoCase( ctx, InlineKeys[i].kname ) ) { |
---|
853 | key = InlineKeys[i].key; |
---|
854 | c = &InlineCode[key]; |
---|
855 | type = InlineKeys[i].type; |
---|
856 | break; |
---|
857 | } |
---|
858 | if( !c ) { |
---|
859 | printf( "\n'%s': Unknown inline option (ignored)", ctx ); |
---|
860 | return; |
---|
861 | } |
---|
862 | |
---|
863 | /* mz_rs_20050607+ */ |
---|
864 | if (c->code) |
---|
865 | totallength = strlen( c->code )+strlen( s ); |
---|
866 | else |
---|
867 | totallength = strlen( s ); |
---|
868 | if (totallength>MAX_INLINE) |
---|
869 | ScanError("\nInline code for %s is too long (%d>%d).\nIncrease MAX_INLINE in scan.h and recompile kpp!", |
---|
870 | ctx, totallength, MAX_INLINE); |
---|
871 | /* mz_rs_20050607- */ |
---|
872 | |
---|
873 | switch( type ) { |
---|
874 | case APPEND: c->code = AppendString( c->code, s, &c->maxlen, MAX_INLINE ); |
---|
875 | break; |
---|
876 | case REPLACE: c->code = ReplaceString( c->code, s, &c->maxlen, MAX_INLINE ); |
---|
877 | break; |
---|
878 | } |
---|
879 | } |
---|
880 | |
---|
881 | int ParseEquationFile( char * filename ) |
---|
882 | { |
---|
883 | int i,j; |
---|
884 | int code; |
---|
885 | |
---|
886 | for( i = 0; i < MAX_SPECIES; i++ ) { |
---|
887 | ReverseCode[i] = NO_CODE; |
---|
888 | Reactive[i] = 0; |
---|
889 | } |
---|
890 | for( i = 0; i < MAX_SPECIES; i++ ) { |
---|
891 | for( j = 0; j < MAX_EQN; j++ ) { |
---|
892 | Stoich_Left[i][j] = 0; |
---|
893 | Stoich[i][j] = 0; |
---|
894 | Stoich_Right[i][j] = 0; |
---|
895 | } |
---|
896 | } |
---|
897 | for( i = 0; i < MAX_SPECIES; i++ ) { |
---|
898 | SpeciesTable[ i ].nratoms = 0; |
---|
899 | } |
---|
900 | |
---|
901 | for( i = 0; i < INLINE_OPT; i++ ) { |
---|
902 | InlineCode[i].code = NULL; |
---|
903 | InlineCode[i].maxlen = 0; |
---|
904 | } |
---|
905 | |
---|
906 | EqnNr = 0; |
---|
907 | SpcNr = 0; |
---|
908 | |
---|
909 | DeclareAtom( "CANCEL" ); |
---|
910 | SetAtomType( "CANCEL", CANCEL_CHECK ); |
---|
911 | DeclareAtom( "IGNORE" ); |
---|
912 | SetAtomType( "IGNORE", NO_CHECK ); |
---|
913 | DeclareSpecies( DUMMY_SPC, "???" ); |
---|
914 | StoreSpecies( MAX_SPECIES-1, DUMMY_SPC, "HV" ); |
---|
915 | AddAtom( "CANCEL", "1" ); |
---|
916 | StoreSpecies( MAX_SPECIES-2, DUMMY_SPC, "PROD" ); |
---|
917 | |
---|
918 | code = Parser( filename ); |
---|
919 | |
---|
920 | return code; |
---|
921 | } |
---|
922 | |
---|