source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/util/mex.f @ 4127

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

Merge of branch palm4u into trunk

File size: 11.8 KB
Line 
1/* #include <math.h> */
2#include "mex.h"
3#include <stdarg.h>
4
5#define MAX_BUF 200
6
7void Usage()
8{
9  mexPrintf("                                                                    \n"
10            "To get this help message use:                KPP_ROOT ?            \n"
11            "                                                                    \n"
12            "To initialize default values use:            KPP_ROOT              \n"
13            "   (type who to see the variables created)                          \n"
14            "                                                                    \n"
15            "To integrate the model use:                                         \n"
16            "          [ c, m, f ] = KPP_ROOT( t, c0, k, p, fn, tfn );          \n"
17            "                                                                    \n"
18            "  input :                                                           \n"
19            "         t   - Time vector, contains the time at which results      \n" 
20            "               should be reported;                                  \n"
21            "         c0  - Vector with the initial concentrations for all       \n"
22            "               species;                                             \n"
23            "         k   - Vector with all rate constants;                      \n"
24            "         p   - Vector of parameters for the integration;            \n"
25            "                 p(1) holds the relative tolerance                  \n"
26            "                 p(2) holds the absolute tolerance                  \n"
27            "                 p(3) holds the minimum step size allowed           \n"
28            "                 p(4) holds the maximum step size allowed           \n"
29            "               If any of the above is zero the default value is     \n"
30            "               used;                                                \n"
31            "         fn  - (optional) Name of a matlab function to be called    \n"
32            "               to update the values of k's and concentrations       \n"                   
33            "               If not present no update is performed.               \n"
34            "                                                                    \n"
35            "         tfn - (optional) Time at which the fn function should      \n"
36            "               be called. If missing <t> is assumed.                \n"
37            "                                                                    \n"
38            "  output:                                                           \n"
39            "         c   - Matrix of concentrations of all species vs. time;    \n"
40            "         m   - (optional) Mass conservation of all atoms vs. time;  \n"
41            "         f   - (optional) Matrix of fluxes of all species vs. time; \n"
42            "                                                                    \n"
43           );
44}
45
46int giveusage;
47
48void F9Error( char *fmt, ... )
49{
50va_list args;
51char buf[ MAX_BUF ];
52char errmsg[ MAX_BUF ];
53
54  va_start( args, fmt );
55  vsprintf( buf, fmt, args );
56  va_end( args );
57 
58  if( giveusage ) Usage();
59 
60  mexPrintf("Error: %s\n", buf); 
61  mexErrMsgTxt( 0 );
62}
63
64char allvars[1000];
65
66int CreateVar(char *name, KPP_REAL val)
67{
68mxArray *GA;
69double  *pga;
70
71  sprintf(allvars, "%s %s",allvars, name);
72  GA = mxCreateDoubleMatrix(1,1,mxREAL);
73  pga = mxGetPr(GA);
74  *pga = (double)val;
75  mxSetName(GA,name);
76  return mexPutArray(GA,"global");
77}
78
79int CreateVec(char *name, int len, KPP_REAL *val)
80{
81mxArray *GA;
82double  *pga;
83int i;
84
85  sprintf(allvars, "%s %s",allvars, name);
86  GA = mxCreateDoubleMatrix(1,len,mxREAL);
87  pga = mxGetPr(GA);
88  if( sizeof(KPP_REAL) == sizeof(double) ) {
89    memmove( pga, val, len*sizeof(double) );
90  } else {
91    for( i = 0; i < len; i++ ) pga[i] = (double)val[i];
92  }   
93  mxSetName(GA,name);
94  return mexPutArray(GA,"global");
95}
96
97int CreateStrVec(char *name, int len, char **val)
98{
99mxArray *GA;
100
101  sprintf(allvars, "%s %s",allvars, name);
102  GA = mxCreateCharMatrixFromStrings( len, (const char **)val );
103  mxSetName(GA,name);
104  return mexPutArray(GA,"global");
105}
106
107int CreateStr(char *name, char *val)
108{
109mxArray *GA;
110
111  sprintf(allvars, "%s %s",allvars, name);
112  GA = mxCreateString( val );
113  mxSetName(GA,name);
114  return mexPutArray(GA,"global");
115}
116
117#define   T_PRM           prhs[0]
118#define  C0_PRM           prhs[1]
119#define   K_PRM           prhs[2]
120#define   P_PRM           prhs[3]
121#define  FN_PRM           prhs[4]
122#define TFN_PRM           prhs[5]
123
124#define   C_PRM           plhs[0]
125#define   M_PRM           plhs[1]
126#define   F_PRM           plhs[2]
127
128#define  HAS_FN           (nrhs >= 5)
129#define  HAS_TFN          (nrhs >= 6)
130
131#define  HAS_M            (nlhs >= 2)
132#define  HAS_F            (nlhs >= 3)
133
134#define  DBL  (sizeof(KPP_REAL) == sizeof(double))
135
136void mexFunction(
137                 int nlhs,       mxArray *plhs[],
138                 int nrhs, const mxArray *prhs[]
139                 )
140{
141double *  tp; 
142double * c0p;
143double *  kp;
144double *  pp;
145char     fnp[ MAX_BUF ];
146double *tfnp;
147
148double *  cp;
149double *  mp;
150double *  fp;
151double ATOLS;
152double  dval[ NMASS+NSPEC ];
153
154mxArray *Carr, *Karr, *Tarr;
155double  *fcp;
156double  *fkp;
157double  *ftp, *ftp1;
158
159int i,j,m,n,nd,t;
160int nsteps, nspc, nreact, ncb;
161int tcb, CallBack;
162KPP_REAL prm[4];
163
164  giveusage = 1;
165
166  if(nrhs == 0) {
167
168    InitVal();
169    Update_RCONST();
170   
171    prm[0] = 1e-4; 
172    prm[1] = 1.0E-18;
173    prm[2] = 0.01;
174    prm[3] = 900;
175
176    sprintf(allvars,"global ");
177   
178    CreateVec("PRM",4, prm);
179   
180    CreateVar("NSPEC",NSPEC);
181    CreateVar("NREACT",NREACT);
182    CreateVar("NMASS",NMASS);
183 
184    CreateVec("C0", NSPEC, C);
185    CreateVec("K0", NREACT, RCONST);
186
187    for( i = 0; i < NLOOKAT; i++ ) 
188      CreateVar( SLOOKAT[i], (double)(i+1) );
189   
190    for( i = 0; i < NMASS; i++ ) 
191      CreateVar( SMASS[i], (double)(i+1) );
192
193    CreateStrVec("SSPEC", NSPEC, SLOOKAT);
194    CreateStrVec("SMASS", NMASS, SMASS);
195    CreateStrVec("SEQN", NREACT, SEQN);
196
197    CreateStr("GLOBALCMD", allvars);
198 
199    mexEvalString(allvars);
200/*
201    mexPrintf("The KPP_ROOT model parameters were sucessfully initialized.\n");   
202*/
203    return;
204  }
205
206  if( nrhs < 4 ) 
207    F9Error("First 4 parameters are REQUIRED only %d received.", nrhs);
208  if( nlhs < 1 ) 
209    F9Error("At least one output parameter REQUIRED.");
210
211  if(! mxIsDouble(T_PRM))   F9Error("<t> must be of type double.");
212  if(! mxIsDouble(C0_PRM))  F9Error("<c0> must be of type double.");
213  if(! mxIsDouble(K_PRM))   F9Error("<k> must be of type double.");
214  if(! mxIsDouble(P_PRM))   F9Error("<p> must be of type double.");
215  if((nrhs > 4) && (! mxIsChar(FN_PRM))) F9Error("<fn> must be of type char.");     
216  if((nrhs > 5) && (! mxIsDouble(TFN_PRM))) F9Error("<tfn> must be of type double.");     
217     
218  nd = mxGetNumberOfDimensions( T_PRM );
219  m  =                  mxGetM( T_PRM ); 
220  n  =                  mxGetN( T_PRM );   
221  if( !( (nd == 2) && ((m == 1) || (n == 1)) ) ) F9Error("<t> must be a column vector.");
222  nsteps = (m == 1) ? n : m;
223  tp = mxGetPr( T_PRM ); 
224
225  nd = mxGetNumberOfDimensions( C0_PRM );
226  m  =                  mxGetM( C0_PRM ); 
227  n  =                  mxGetN( C0_PRM );   
228  if( !( (nd == 2) && ((m == 1) || (n == 1)) ) ) F9Error("<c0> must be a column vector.");
229  nspc = (m == 1) ? n : m;
230  c0p = mxGetPr( C0_PRM ); 
231
232  nd = mxGetNumberOfDimensions( K_PRM );
233  m  =                  mxGetM( K_PRM ); 
234  n  =                  mxGetN( K_PRM );   
235  if( !( (nd == 2) && ((m == 1) || (n == 1)) ) ) F9Error("<k> must be a column vector.");
236  nreact = (m == 1) ? n : m;
237  kp = mxGetPr( K_PRM ); 
238
239  nd = mxGetNumberOfDimensions( P_PRM );
240  m  =                  mxGetM( P_PRM ); 
241  n  =                  mxGetN( P_PRM );   
242  if( !( (nd == 2) && ((m == 1) || (n == 1)) && (n*m == 4) ) )
243    F9Error("<p> must be a column vectorof length 4.");
244  pp = mxGetPr( P_PRM ); 
245
246  *fnp = 0;
247  if( HAS_FN ) {
248    nd = mxGetNumberOfDimensions( FN_PRM );
249    m  =                  mxGetM( FN_PRM ); 
250    n  =                  mxGetN( FN_PRM );   
251    if( !( (nd == 2) && ((m == 1) || (n == 1)) ) ) F9Error("<fn> must be a character string.");
252    if( mxGetString( FN_PRM, fnp, MAX_BUF ) ) 
253      F9Error("Can not read function mane (too long?)"); 
254
255    Carr = mxCreateDoubleMatrix(1,NSPEC,mxREAL);
256    fcp = mxGetPr(Carr);
257    mxSetName(Carr,"C");
258    mexPutArray(Carr,"base");
259   
260    Karr = mxCreateDoubleMatrix(1,NREACT,mxREAL);
261    fkp = mxGetPr(Karr);
262    mxSetName(Karr,"K");
263    mexPutArray(Karr,"base");
264   
265    Tarr = mxCreateDoubleMatrix(1,1,mxREAL);
266    ftp = mxGetPr(Tarr);
267    mxSetName(Tarr,"T");
268    mexPutArray(Tarr,"base");
269  }
270
271  tfnp = 0; ncb = 0;
272  if( HAS_TFN ) {
273    nd = mxGetNumberOfDimensions( TFN_PRM );
274    m  =                  mxGetM( TFN_PRM ); 
275    n  =                  mxGetN( TFN_PRM );   
276    if( !( (nd == 2) && ((m == 1) || (n == 1)) ) ) F9Error("<tfn> must be a column vector.");
277    ncb = (m == 1) ? n : m; 
278    tfnp = mxGetPr( TFN_PRM );
279  }
280
281  giveusage = 0;
282
283  if( !((nspc == NSPEC) && (nreact == NREACT)) ) {
284    F9Error("Size of parameters do not match the model:\n\n"
285            "  Number of species was %d and should be %d;\n"             
286            "  Number of rections (rate constants) was %d and should be %d;\n",
287            nspc, NSPEC, nreact, NREACT);   
288  }
289
290  if( DBL ) { memmove( C, c0p, sizeof(double)*NSPEC );
291              memmove( RCONST, kp, sizeof(double)*NREACT ); }
292       else { for( i = 0; i < NSPEC; i++ ) C[i] = (KPP_REAL)c0p[i];
293              for( i = 0; i < NREACT; i++ ) RCONST[i] = (KPP_REAL)kp[i]; }
294 
295  RTOLS = 1e-4;
296  ATOLS = 1e-18;
297  STEPMIN = 0.01;
298  STEPMAX = 900.0;
299         
300  if( pp[0] ) RTOLS = pp[0];
301  if( pp[1] ) ATOLS = pp[1];
302  if( pp[2] ) STEPMIN = pp[2];
303  if( pp[3] ) STEPMAX = pp[3];
304 
305  for( i = 0; i < NVAR; i++ ) {
306    RTOL[i] = RTOLS;
307    ATOL[i] = ATOLS;
308  }
309 
310  C_PRM = mxCreateDoubleMatrix(NSPEC,nsteps,mxREAL); 
311  cp = mxGetPr(C_PRM);
312 
313  if( HAS_M ) {
314    M_PRM = mxCreateDoubleMatrix(NMASS,nsteps,mxREAL);
315    mp = mxGetPr(M_PRM);
316  }
317 
318  if( HAS_F ) {
319    F_PRM = mxCreateDoubleMatrix(NSPEC,nsteps,mxREAL);
320    fp = mxGetPr(F_PRM);
321  }
322
323  tcb = 0;
324
325  for( t = 0; t < nsteps; t++ ) {             
326    if( t ) {
327      TIME = tp[t-1];
328
329      CallBack = 0;
330      if( HAS_TFN ) {
331        if( tcb < ncb )
332          if( tfnp[tcb] <= TIME ) { CallBack = 1; tcb++; }
333      } else {
334        CallBack = HAS_FN;
335      } 
336
337      if( CallBack ) {
338        if( DBL ) { memmove( fcp, C, sizeof(double)*NSPEC ); 
339                    memmove( fkp, RCONST, sizeof(double)*NREACT ); }
340             else { for( i = 0; i < NSPEC; i++ ) fcp[i] = (double)C[i]; 
341                    for( i = 0; i < NREACT; i++ ) fkp[i] = (double)RCONST[i]; }
342        *ftp = TIME; 
343
344        mexPutArray(Carr,"base");
345        mexPutArray(Karr,"base");
346        mexPutArray(Tarr,"base");
347
348        mexCallMATLAB( 0, 0, 0, 0, fnp );
349
350        mxDestroyArray(Carr); Carr = mexGetArray("C","base"); fcp = mxGetPr(Carr);
351        mxDestroyArray(Karr); Karr = mexGetArray("K","base"); fkp = mxGetPr(Karr);
352        mxDestroyArray(Tarr); Tarr = mexGetArray("T","base"); ftp = mxGetPr(Tarr);
353
354        if( DBL ) { memmove( C, fcp, sizeof(double)*NSPEC ); 
355                    memmove( RCONST, fkp, sizeof(double)*NREACT ); }
356             else { for( i = 0; i < NSPEC; i++ ) C[i] = (KPP_REAL)fcp[i]; 
357                    for( i = 0; i < NREACT; i++ ) RCONST[i] = (KPP_REAL)fkp[i]; }
358 
359      }
360
361      INTEGRATE( tp[t-1], tp[t] );
362    }
363    if( DBL ) { memmove( cp, C, sizeof(double)*NSPEC ); cp += NSPEC; } 
364         else { for( i = 0; i < NSPEC; i++ ) *cp++ = (double)C[i]; }
365    if( HAS_M ) {
366      if( DBL ) { GetMass( mp ); mp += NMASS; } 
367           else { GetMass( dval );
368                  for( i = 0; i < NMASS; i++ ) *mp++ = (double)dval[i]; }
369    }
370    if( HAS_F ) {
371      if( DBL ) { FLUX( fp ); fp += NSPEC; }
372           else { FLUX( dval );
373                  for( i = 0; i < NSPEC; i++ ) *fp++ = (double)dval[i]; }           
374    }
375  }
376}
377
Note: See TracBrowser for help on using the repository browser.