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

Last change on this file since 4843 was 4843, checked in by raasch, 16 months ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

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) -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 <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.