source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/src/kpp.c @ 3225

Last change on this file since 3225 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

File size: 14.2 KB
Line 
1/******************************************************************************
2
3  KPP - The Kinetic PreProcessor
4        Builds simulation code for chemical kinetic systems
5
6  Copyright (C) 1995-1996 Valeriu Damian and Adrian Sandu
7  Copyright (C) 1997-2005 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 <stdlib.h>
34#include <string.h>
35#include "gdata.h"
36#include "scan.h"
37
38char *eqFileName;
39char *rootFileName = "ff";
40char Home[ MAX_PATH ] = ""; 
41
42short int linStru[ MAX_SPECIES ];
43short int colStru[ MAX_SPECIES ];
44short int bestStru[ MAX_SPECIES ];
45short int *Stru;
46
47enum stru_criteria { UNSORT, LINSORT, COLSORT, BESTSORT };
48
49void EqCopy( EQ_VECT e1, EQ_VECT e2 )
50{
51int i;
52
53  for( i = 0; i < EqnNr; i++ ) e2[i] = e1[i];
54}
55
56int NoSort( const void *p1, const void *p2 )
57{
58  return -1;
59}
60
61int CodeCmp( const void *p1, const void *p2 )
62{
63CODE *c1, *c2;
64
65  c1 = (CODE*)p1; 
66  c2 = (CODE*)p2; 
67 
68  if ( *c1 < *c2 ) return -1;
69  if ( *c1 > *c2 ) return 1;
70  return 0;
71}
72
73int CodeRCmp( const void *p1, const void *p2 )
74{
75int rc1, rc2;
76CODE *c1, *c2;
77
78  c1 = (CODE*)p1; 
79  c2 = (CODE*)p2; 
80
81  rc1 = Reactive[ ReverseCode[ *c1 ] ]; 
82  rc2 = Reactive[ ReverseCode[ *c2 ] ];
83  if ( rc1 > rc2 ) return -1; 
84  if ( rc1 < rc2 ) return 1;
85  if ( *c1 < *c2 ) return -1;
86  if ( *c1 > *c2 ) return 1;
87  return 0;
88}
89
90int CodeSSCmp( const void *p1, const void *p2 )
91{
92  return -CodeRCmp(p1,p2);
93}
94
95int CodeSCmp( const void *p1, const void *p2 )
96{
97CODE *c1, *c2;
98short int sc1, sc2;
99
100  c1 = (CODE*)p1; 
101  c2 = (CODE*)p2; 
102 
103  sc1 = Stru[ ReverseCode[ *c1 ] ];
104  sc2 = Stru[ ReverseCode[ *c2 ] ];
105 
106  if ( sc1 > sc2 ) return 1; 
107  if ( sc1 < sc2 ) return -1;
108  if ( *c1 < *c2 ) return 1;
109  if ( *c1 > *c2 ) return -1;
110  return 0;
111}
112
113void UpdateStructJ()
114{
115int i,j,k;
116
117  for ( i=0; i<VarNr; i++ )
118    for ( j=0; j<VarNr; j++ )
119      structJ[i][j]=(i==j)?1:0;
120             
121  for (i = 0; i < VarNr; i++)
122    for (j = 0; j < VarNr; j++)
123      for (k = 0; k < EqnNr; k++)
124        if ( Stoich[i][k]*((Stoich_Left[j][k])?1:0) != 0.0 )
125          structJ[i][j]=1;
126
127  for ( i=0; i<VarNr; i++ )
128    for ( j=0; j<VarNr; j++ )
129      LUstructJ[i][j]=structJ[i][j];
130             
131}
132
133int ComputeLUStructJ()
134{
135int i,j,k;
136int nu,nl;
137
138  for (j = 0; j < VarNr-1; j++) {
139    for (i = j+1; i < VarNr; i++) {
140      if( LUstructJ[i][j] ) {
141        for (k = j; k < VarNr; k++) 
142         /*   LUstructJ[i][k] += LUstructJ[j][k]; */
143           if ( LUstructJ[j][k] != 0 )
144              LUstructJ[i][k] = 1; 
145      }   
146    }
147  } 
148
149 
150  nu = 0; nl = 0;
151  for (i = 0; i < VarNr; i++) 
152    for (j = 0; j < VarNr; j++)
153      if( LUstructJ[i][j] ) {
154        if(i > j) nl++;
155        if(i <= j) nu++;     
156      }         
157
158  return nu+nl;
159}
160
161int LUnonZero()
162{
163CODE v[MAX_SPECIES];
164CODE *var;
165int i,j,k;
166int nu,nl;
167
168  var = v;
169  if( Stru != bestStru ) {
170    for( i=0; i<VarNr; i++ )
171      var[i] = Code[i];
172    qsort( (void*)var, VarNr, sizeof(CODE), CodeSCmp );
173  } else {
174    var = bestStru;
175  }
176         
177  for (i = 0; i < VarNr; i++)
178    for (j = 0; j < VarNr; j++)
179      LUstructJ[i][j] = structJ[ ReverseCode[var[i]] ][ ReverseCode[var[j]] ]; 
180
181  return ComputeLUStructJ();
182}
183
184void LinColSparsity()
185{
186int i,j,k;
187int nlin, ncol;
188FILE * fff;
189
190  for ( i=0; i<VarNr; i++ )
191    for ( j=0; j<VarNr; j++ )
192      structJ[i][j]=(i==j)?1:0;
193             
194  for (i = 0; i < VarNr; i++)
195    for (j = 0; j < VarNr; j++)
196      for (k = 0; k < EqnNr; k++)
197        if ( Stoich[i][k]*((Stoich_Left[j][k])?1:0) != 0.0 )
198          structJ[i][j]=1;
199
200  for ( i=0; i<VarNr; i++ ) {
201    linStru[i] = 0;
202    for (j = 0; j < VarNr; j++)
203      linStru[i] += structJ[i][j];
204  }
205
206  for ( i=0; i<VarNr; i++ ) {
207    colStru[i] = 0;
208    for (j = 0; j < VarNr; j++)
209      colStru[i] += structJ[j][i];
210    colStru[i] *= linStru[i]; 
211  }
212
213  Stru = linStru;
214  nlin = LUnonZero();
215  Stru = colStru;
216  ncol = LUnonZero();
217  if( nlin <= ncol ) { 
218    Stru = linStru;
219    LUnonZero();
220  } 
221}
222
223void BestSparsity()
224{
225int i,j,k;
226int cnz, lnz;
227int best, crt;
228int best_i;
229int tmp;
230int s;
231
232  UpdateStructJ();
233     
234  for ( i=0; i<VarNr; i++ )
235    bestStru[i] = Code[i];
236
237  for ( s=0; s<VarNr-1; s++ ) {
238    best = MAX_SPECIES*MAX_SPECIES; best_i = 0;
239    for ( i=s; i<VarNr; i++ ) {
240      cnz = 0;lnz = 0;
241      for (j = s; j < VarNr; j++) {
242        cnz += (LUstructJ[i][j]?1:0);
243        lnz += (LUstructJ[j][i]?1:0);
244      }
245      crt = (cnz-1)*(lnz-1);
246      if( crt < best ) {
247        best = crt;
248        best_i = i;
249      }   
250    }
251    for ( i=0; i<VarNr; i++ ) {
252      tmp = LUstructJ[s][i];
253      LUstructJ[s][i] = LUstructJ[best_i][i];
254      LUstructJ[best_i][i] = tmp;     
255    }   
256    for ( i=0; i<VarNr; i++ ) {
257      tmp = LUstructJ[i][s];
258      LUstructJ[i][s] = LUstructJ[i][best_i];
259      LUstructJ[i][best_i] = tmp;     
260    }   
261    tmp = bestStru[s];
262    bestStru[s] = bestStru[best_i];
263    bestStru[best_i] = tmp;
264   
265    for (i = s+1; i < VarNr; i++) {
266      if( LUstructJ[i][s] ) {
267        for (k = s; k < VarNr; k++) 
268          LUstructJ[i][k] += LUstructJ[s][k];
269      }   
270    }
271  }
272
273  Stru = bestStru;
274}
275
276void ReorderSpecies( int criteria )
277{
278CODE *var;
279CODE *fix;
280CODE *dummy;
281EQ_VECT *tmpStoich_Left;
282EQ_VECT *tmpStoich_Right;
283EQ_VECT *tmpStoich;
284CODE *tmpCode;
285CODE *tmpReact;
286int i, k;
287int new;
288int (*cmpVar)(const void *, const void *);
289int (*cmpFix)(const void *, const void *);
290int dummyNr;
291
292  cmpVar = CodeCmp;
293  cmpFix = CodeCmp;
294
295  switch( criteria ) {
296    case UNSORT:   cmpVar = useJacobian ? CodeRCmp : CodeCmp;
297                   break; 
298    case LINSORT:  cmpVar = useJacobian ? CodeSCmp : CodeCmp;
299                   Stru = linStru; 
300                   break; 
301    case COLSORT:  cmpVar = useJacobian ? CodeSCmp : CodeCmp;
302                   Stru = colStru;
303                   break; 
304    case BESTSORT: cmpVar = useJacobian ? NoSort : CodeCmp;
305                   break; 
306  }
307
308  VarNr = 0;
309  VarActiveNr = 0;
310  FixNr = 0;
311  dummyNr = 0;
312
313  var = (CODE*)malloc( SpcNr * sizeof(CODE) );
314  fix = (CODE*)malloc( SpcNr * sizeof(CODE) );
315  dummy = (CODE*)malloc( 5 * sizeof(CODE) );
316  tmpStoich_Left = (EQ_VECT*)malloc( SpcNr * sizeof(EQ_VECT) );
317  tmpStoich_Right = (EQ_VECT*)malloc( SpcNr * sizeof(EQ_VECT) );
318  tmpStoich = (EQ_VECT*)malloc( SpcNr * sizeof(EQ_VECT) );
319  tmpCode = (CODE*)malloc( SpcNr * sizeof(CODE) );
320  tmpReact = (CODE*)malloc( SpcNr * sizeof(CODE) );
321
322  for( i = 0; i < SpcNr; i++ ) {
323    switch( SpeciesTable[ Code[i] ].type ) {
324      case VAR_SPC:  var[ VarNr++ ] = Code[ i ];
325                     break;
326      case FIX_SPC:  fix[ FixNr++ ] = Code[ i ]; 
327                     break;
328      case DUMMY_SPC:dummy[ dummyNr++ ] = Code[ i ]; 
329                     break;
330    }
331  }
332
333  if( Stru != bestStru ) {
334    qsort( (void*)var, VarNr, sizeof(CODE), cmpVar );
335  } else {
336    for( i = 0; i < SpcNr; i++ )
337      var[i] = bestStru[i];
338  } 
339  qsort( (void*)fix, FixNr, sizeof(CODE), cmpFix );
340
341  for( i = 0; i < SpcNr; i++ ) {
342    EqCopy( Stoich_Left[i], tmpStoich_Left[i] );
343    EqCopy( Stoich_Right[i], tmpStoich_Right[i] );
344    EqCopy( Stoich[i], tmpStoich[i] );
345    tmpCode[i] = Code[i];
346    tmpReact[i] = Reactive[i];
347  }
348
349  SpcNr -= dummyNr;
350  dummyNr = 0;
351
352  k = 0;
353  for( i = 0; i < VarNr; i++ ) {
354    new = ReverseCode[ var[i] ];
355    EqCopy( tmpStoich_Left[ new ], Stoich_Left[ k ] );
356    EqCopy( tmpStoich_Right[ new ], Stoich_Right[ k ] );
357    EqCopy( tmpStoich[ new ], Stoich[ k ] );
358    Code[ k ] = tmpCode[ new ];
359    Reactive[ k ] = tmpReact[ new ];
360    if( Reactive[ k ] ) VarActiveNr++; 
361    k++;
362  }
363  for( i = 0; i < FixNr; i++ ) {
364    new = ReverseCode[ fix[i] ];
365    EqCopy( tmpStoich_Left[ new ], Stoich_Left[ k ] );
366    EqCopy( tmpStoich_Right[ new ], Stoich_Right[ k ] );
367    EqCopy( tmpStoich[ new ], Stoich[ k ] );
368    Code[ k ] = tmpCode[ new ];
369    Reactive[ k ] = tmpReact[ new ];
370    k++;
371  }
372  for( i = 0; i < dummyNr; i++ ) {
373    new = ReverseCode[ dummy[i] ];
374    EqCopy( tmpStoich_Left[ new ], Stoich_Left[ k ] );
375    EqCopy( tmpStoich_Right[ new ], Stoich_Right[ k ] );
376    EqCopy( tmpStoich[ new ], Stoich[ k ] );
377    Code[ k ] = tmpCode[ new ];
378    Reactive[ k ] = tmpReact[ new ];
379    k++;
380  }
381
382
383  for( i = 0; i < SpcNr+dummyNr; i++ )
384    ReverseCode[ Code[i] ] = i;
385
386  free( tmpReact );
387  free( tmpCode );
388  free( tmpStoich );
389  free( tmpStoich_Right );
390  free( tmpStoich_Left );
391  free( dummy );
392  free( fix );
393  free( var );   
394
395  fflush(stdout);
396}
397
398/* Allocate Internal Arrays */
399void  AllocInternalArrays( void )
400{
401int i;
402
403if ( (Stoich_Left =(float**)calloc(MAX_SPECIES,sizeof(float*)))==NULL ) 
404    FatalError(-30,"Cannot allocate Stoich_Left.\n");
405
406for (i=0; i<MAX_SPECIES; i++)   
407    if ( (Stoich_Left[i] = (float*)calloc(MAX_EQN,sizeof(float)))==NULL ) {
408        FatalError(-30,"Cannot allocate Stoich_Left[%d]",i,MAX_SPECIES);
409    }
410
411if ( (Stoich_Right = (float**)calloc(MAX_SPECIES,sizeof(float*)))==NULL ) 
412    FatalError(-30,"Cannot allocate Stoich_Right.\n");
413
414for (i=0; i<MAX_SPECIES; i++)   
415    if ( (Stoich_Right[i] = (float*)calloc(MAX_EQN,sizeof(float)))==NULL ) {
416        FatalError(-30,"Cannot allocate Stoich_Right[%d].",i);
417    }
418
419if ( (Stoich = (float**)calloc(MAX_SPECIES,sizeof(float*)))==NULL ) 
420    FatalError(-30,"Cannot allocate Stoich.\n");
421
422for (i=0; i<MAX_SPECIES; i++)   
423    if ( (Stoich[i] = (float*)calloc(MAX_EQN,sizeof(float)))==NULL ) {
424        FatalError(-30,"Cannot allocate Stoich[%d].",i);
425    }
426
427}
428
429
430/* Allocate Structure Arrays */
431void  AllocStructArrays( void )
432{
433int i;
434
435
436if ( (structB = (int**)calloc(EqnNr,sizeof(int*)))==NULL ) 
437    FatalError(-30, "Cannot allocate structB.");
438
439for (i=0; i<EqnNr; i++)   
440    if ( (structB[i] =(int*) calloc(SpcNr,sizeof(int)))==NULL )
441        FatalError(-30, "Cannot allocate structB[%d].\n",i);
442   
443if ( (structJ = (int**)calloc(SpcNr,sizeof(int*)))==NULL ) 
444    FatalError(-30, "Cannot allocate structJ.");
445
446for (i=0; i<SpcNr; i++)   
447    if ( (structJ[i] =(int*) calloc(SpcNr,sizeof(int)))==NULL ) 
448        FatalError(-30, "Cannot allocate structJ[%d].\n",i);
449   
450if ( (LUstructJ = (int**)calloc(SpcNr,sizeof(int*)))==NULL ) 
451    FatalError(-30, "Cannot allocate LUstructJ.");
452
453for (i=0; i<SpcNr; i++)   
454    if ( (LUstructJ[i] = (int*)calloc(SpcNr,sizeof(int)))==NULL ) 
455        FatalError(-30, "Cannot allocate LUstructJ[%d].\n",i);
456
457}
458
459/*******************************************************************/                   
460int Postprocess( char * root )
461{
462char buf[ 200 ];
463char cmd[500];
464char cmdexe[500];
465static char tmpfile[] = "kppfile.tmp";
466FILE * fp;
467
468  if ( useLang == MATLAB_LANG ) {
469 /*  Add rate function definitions as internal functions to the Update_RCONST file*/
470    sprintf( buf, "cat %s_Update_RCONST.m %s_Rates.m > tmp; mv tmp %s_Update_RCONST.m;", 
471        root, root, root ); 
472    system( buf );       
473  }
474
475/*    Postprocessing to replace parameter names by values in the declarations
476  strcpy( cmd, "sed " );
477  sprintf( cmd, "%s -e 's/(NVAR)/(%d)/g'", cmd, VarNr ); 
478  sprintf( cmd, "%s -e 's/(NFIX)/(%d)/g'", cmd, FixNr ); 
479  sprintf( cmd, "%s -e 's/(NSPEC)/(%d)/g'", cmd,SpcNr ); 
480  sprintf( cmd, "%s -e 's/(NREACT)/(%d)/g'", cmd, EqnNr ); 
481  sprintf( cmd, "%s -e 's/(NONZERO)/(%d)/g'", cmd, Jac_NZ ); 
482  sprintf( cmd, "%s -e 's/(LU_NONZERO)/(%d)/g'", cmd, LU_Jac_NZ ); 
483  sprintf( cmd, "%s -e 's/(NHESS)/(%)/g'", cmd, Hess_NZ ); 
484   
485  sprintf( buf, "%s_Function", rootFileName ); 
486  switch( useLang ) {
487    case F77_LANG: sprintf( buf, "%s.f", buf );
488                 break;
489    case F90_LANG: sprintf( buf, "%s.f90", buf );
490                 break;
491    case C_LANG: sprintf( buf, "%s.c", buf );
492                 break;
493    case MATLAB_LANG: sprintf( buf, "%s.m", buf );
494                 break;
495    default: printf("\n Language '%d' not implemented!\n",useLang);
496                 exit(1);
497  }
498  sprintf( cmdexe, "%s %s > %s; mv %s %s;", cmd, buf, tmpfile, tmpfile, buf ); 
499  printf("\n\nCMDEXE='%s'\n",cmdexe);
500  system( cmdexe );
501*/
502} 
503 
504/*******************************************************************/                   
505int main( int argc, char * argv[] )
506{
507int status;
508char name[ 200 ];
509char *p;
510int i,j;
511 
512  AllocInternalArrays();
513
514  p = getenv("KPP_HOME");
515  if( p ) strcpy( Home, p );
516
517  switch( argc ) {
518    case 3: eqFileName = argv[1];
519            rootFileName = argv[2];
520            break;
521    case 2: eqFileName = argv[1];
522            strcpy( name, eqFileName );
523            p = name + strlen(name);
524            while( p > name ) {
525              if( *p == '.') { 
526                *p = '\0';
527                break;
528              }
529              p--; 
530            } 
531            rootFileName = name;
532            break;
533    default: FatalError(1,"\nUsage :"
534                          "\n        kpp <equations file> [output file]\n");
535  }
536
537  printf("\nThis is KPP-%s.\n", KPP_VERSION);
538
539  printf("\nKPP is parsing the equation file.");
540  status = ParseEquationFile( argv[1] );
541
542  if( status ) FatalError(2,"%d errors and %d warnings encountered.", 
543                           nError, nWarning ); 
544  /* Allocate some internal data structures */
545  AllocStructArrays();
546
547  printf("\nKPP is computing Jacobian sparsity structure.");
548  ReorderSpecies( UNSORT );
549  if (useReorder==1){
550    BestSparsity(); 
551    ReorderSpecies( BESTSORT );
552    }
553  UpdateStructJ();
554  ComputeLUStructJ();
555
556  if( initNr == -1 ) initNr = VarNr;
557
558
559  printf("\nKPP is starting the code generation.");
560  Generate( rootFileName );
561 
562  printf("\nKPP is starting the code post-processing.");
563  Postprocess( rootFileName );
564 
565  printf("\n\nKPP has succesfully created the model \"%s\".\n\n",rootFileName);
566
567  if( nError ) exit(4);
568  if( nWarning ) exit(5);
569 
570  exit(0);
571}
Note: See TracBrowser for help on using the repository browser.