source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/oldies/ros2.c @ 4321

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

Merge of branch palm4u into trunk

File size: 7.7 KB
Line 
1
2        #define MAX(a,b) ((a) >= (b)) ?(a):(b)
3        #define MIN(b,c) ((b) <  (c)) ?(b):(c) 
4        #define abs(x)   ((x) >=  0 ) ?(x):(-x)
5        #define dabs(y)  (double)abs(y)
6        #define DSQRT(d) (double)pow(d,0.5)
7        #define signum(x)((x) >=  0 ) ?(1):(-1)
8
9
10
11        void (*forfun)(int,KPP_REAL,KPP_REAL [],KPP_REAL []);
12        void (*forjac)(int,KPP_REAL,KPP_REAL [],KPP_REAL []);
13
14
15
16 
17void FUNC_CHEM(int N,KPP_REAL T,KPP_REAL Y[NVAR],KPP_REAL P[NVAR])
18        {
19         KPP_REAL Told;
20         Told = TIME;
21         TIME = T;
22         Update_SUN();
23         Update_RCONST();
24         Fun( Y, FIX, RCONST, P );
25         TIME = Told;
26        }
27 
28void JAC_CHEM(int N,KPP_REAL T,KPP_REAL Y[NVAR],KPP_REAL J[LU_NONZERO])
29        {
30         KPP_REAL Told;
31         Told = TIME;
32         TIME = T;
33         Update_SUN();
34         Update_RCONST();
35         Jac_SP( Y, FIX, RCONST, J );
36         TIME = Told;
37        } 
38                                                                                                                         
39
40
41
42
43       INTEGRATE(KPP_REAL TIN,KPP_REAL TOUT )
44        {
45 
46         /*  TIN - Start Time   */
47         /*  TOUT - End Time     */
48
49
50         int INFO[5];
51         forfun = &FUNC_CHEM;
52         forjac = &JAC_CHEM;   
53         INFO[0] = Autonomous;
54         ROS2(NVAR,TIN,TOUT,STEPMIN,STEPMAX,STEPMIN,VAR,ATOL
55         ,RTOL,INFO,forfun,forjac);
56 
57        }
58 
59 
60int ROS2(int N,KPP_REAL T, KPP_REAL Tnext,KPP_REAL Hmin,KPP_REAL Hmax,
61   KPP_REAL Hstart,KPP_REAL y[NVAR],KPP_REAL AbsTol[NVAR],KPP_REAL RelTol[NVAR],
62   int INFO[5],void (*forfun)(int,KPP_REAL,KPP_REAL [],KPP_REAL []),
63   void (*forjac)(int,KPP_REAL,KPP_REAL [],KPP_REAL []) )
64   {
65     
66                                                                               
67/*
68     All the arguments aggree with the KPP syntax.
69
70  INPUT ARGUMENTS:
71     y = Vector of (NVAR) concentrations, contains the
72         initial values on input
73     [T, Tnext] = the integration interval
74     Hmin, Hmax = lower and upper bounds for the selected step-size.
75          Note that for Step = Hmin the current computed
76          solution is unconditionally accepted by the error
77          control mechanism.
78     AbsTol, RelTol = (NVAR) dimensional vectors of
79          componentwise absolute and relative tolerances.
80     FUNC_CHEM = name of routine of derivatives. KPP syntax.
81          See the header below.
82     JAC_CHEM = name of routine that computes the Jacobian, in
83          sparse format. KPP syntax. See the header below.
84     Info(1) = 1  for  autonomous   system
85             = 0  for nonautonomous system
86
87  OUTPUT ARGUMENTS:
88     y = the values of concentrations at Tend.
89     T = equals Tend on output.
90     Info(2) = # of FUNC_CHEM calls.
91     Info(3) = # of JAC_CHEM calls.
92     Info(4) = # of accepted steps.
93     Info(5) = # of rejected steps.
94*/
95 
96      KPP_REAL K1[NVAR], K2[NVAR], K3[NVAR], K4[NVAR];
97      KPP_REAL F1[NVAR], JAC[LU_NONZERO];
98      KPP_REAL ghinv , uround , dround , c43 , tau;
99      KPP_REAL ynew[NVAR];
100      KPP_REAL H, Hold, Tplus;
101      KPP_REAL ERR, factor, facmax;
102      int n,nfcn,njac,Naccept,Nreject,i,j,ier;
103      char IsReject,Autonomous;
104                                                                           
105      KPP_REAL gamma, m1, m2, alpha, beta, delta, theta, g[NVAR], x[NVAR];
106 
107/*     Initialization of counters, etc.       */
108      Autonomous = (INFO[0] == 1);         
109      uround = (double)(1e-15);               
110
111      dround = DSQRT(uround);
112      c43 = (double)(- 8.e0/3.e0);
113      H = MAX( (double)1.e-8, Hmin );
114      Tplus = T;
115      IsReject = 0;
116      Naccept  = 0;
117      Nreject  = 0;
118      nfcn     = 0;
119      njac     = 0;
120      gamma = (double)(1.e0 + 1.e0/DSQRT(2.e0));
121 
122/* === Starting the time loop ===     */
123 while(T < Tnext)
124 {
125 ten :     
126     Tplus = T + H;
127
128     if ( Tplus > Tnext )
129      {
130        H = Tnext - T;
131        Tplus = Tnext;
132      }
133
134     (*forjac)(NVAR, T, y,JAC );               
135
136     njac = njac+1;
137     ghinv = (double)(-1.0e0/(gamma*H));
138     
139     
140     
141     for(j=0;j<NVAR;j++)
142        JAC[LU_DIAG[j]] = JAC[LU_DIAG[j]] + ghinv;
143     
144
145
146     ier = KppDecomp (JAC);
147
148     if ( ier != 0)
149     {
150      if( H > Hmin ) 
151      {
152       H = (double)(5.0e-1*H); 
153       goto ten;           
154      }
155     else
156       printf("IER <> 0 , H = %d", H);
157       
158     }/* main ier if ends*/     
159               
160
161     (*forfun)(NVAR , T, y, F1 ) ;             
162
163
164
165                                                                                                                           
166
167 
168/* ====== NONAUTONOMOUS CASE ===============   */
169       if(Autonomous == 0) 
170       {
171         tau =( dround*MAX ((double)1.0e-6, dabs(T)) *signum(T) );
172         (*forfun)(NVAR, T+tau, y, K2);
173         nfcn=nfcn+1;
174         
175         for(j = 0;j<NVAR;j++)
176             K3[j] = ( K2[j]-F1[j] )/tau;
177         
178 
179/* ----- STAGE 1 (NON-AUTONOMOUS) -----   */
180
181         delta = (double)(gamma*H);
182         
183         for(j = 0;j<NVAR;j++)
184             K1[j] =  F1[j] + delta*K3[j];
185         
186         KppSolve (JAC, K1);
187 
188
189/* ----- STAGE 2  (NON-AUTONOMOUS) -----   */
190         alpha = (double)(- 1.e0/gamma);
191         for(j = 0;j<NVAR;j++)
192             ynew[j] = y[j] + alpha*K1[j];
193         
194 
195         (*forfun)(NVAR, T+H, ynew, F1);
196         nfcn=nfcn+1;
197         beta = (double)(2.e0/(gamma*H));
198         delta = (double)(-gamma*H);
199         for(j = 0;j<NVAR;j++)
200             K2[j] = F1[j] + beta*K1[j] + delta*K3[j];
201         
202         KppSolve (JAC, K2);
203   
204       }/* if for non - autonomous case ends here  */   
205 
206/* ====== AUTONOMOUS CASE ===============    */
207       else
208       {
209 
210/* ----- STAGE 1 (AUTONOMOUS) -----   */
211         for(j = 0;j<NVAR;j++)
212             K1[j] =  F1[j];
213         
214         KppSolve (JAC, K1);
215 
216/* ----- STAGE 2  (AUTONOMOUS) -----  */
217         alpha = (double)(- 1.e0/gamma);
218       
219         for(j = 0;j<NVAR;j++)
220             ynew[j] = y[j] + alpha*K1[j];
221       
222         (*forfun)(NVAR, T+H, ynew, F1);
223
224         nfcn=nfcn+1;
225         beta = (double)(2.e0/(gamma*H));
226         for(j = 0;j < NVAR;j++)
227             K2[j] = F1[j] + beta*K1[j];
228       
229         KppSolve (JAC, K2);
230       
231       }/* else autonomous case ends here */
232       
233 
234                                                                                                           
235/* ---- The Solution ---  */
236       m1 = (double)(-3.e0/(2.e0*gamma));
237       m2 = (double)(-1.e0/(2.e0*gamma));
238       for(j = 0;j<NVAR;j++)
239          ynew[j] = y[j] + m1*K1[j] + m2*K2[j];
240       
241 
242 
243/* ====== Error estimation ========    */
244 
245        ERR=(double)0.e0;
246        theta = (double)(1.e0/(2.e0*gamma));
247       
248        for(i=0;i<NVAR;i++)
249        {
250           g[i] = AbsTol[i] + RelTol[i]*dabs(ynew[i]);
251           ERR = (double)( ERR + pow( ( theta*(K1[i]+K2[i])/g[i] ) , 2.0 ) );
252        }
253         ERR = MAX( uround, DSQRT( ERR/NVAR ) );
254 
255/* ======= Choose the stepsize ==================  */
256 
257 
258        factor = (double)(0.9/pow( ERR,(1.e0/2.e0) ));
259        if (IsReject == 1) 
260            facmax=(double)1.0;
261       
262        else
263            facmax=(double)10.0;
264       
265       
266        factor =(double) MAX( 1.0e-1, MIN(factor,facmax) );
267        Hold = H;
268        H = (double)MIN( Hmax, MAX(Hmin,factor*H) );
269                                                                                               
270/* ======= Rejected/Accepted Step ==================  */
271 
272 
273        if( (ERR>1) && (Hold>Hmin) )
274        {
275         IsReject = 1;
276         Nreject = Nreject + 1;
277        }               
278        else
279        {
280         IsReject = 0;
281         for(i=0;i<NVAR;i++)
282         {
283          y[i] = ynew[i];
284         }
285        T = Tplus;
286        Naccept = Naccept+1;   
287 
288        }/*  else ends here */
289
290 
291 
292 
293 
294/* ======= End of the time loop =================    */
295 
296 }    /*      while loop (T < Tnext) ends here    */ 
297 
298 
299/* ======= Output Information ================      */
300      INFO[1] = nfcn;
301      INFO[2] = njac;
302      INFO[3] = Naccept;
303      INFO[4] = Nreject;
304
305      return 0; 
306       
307  } /* function ros2 ends here   */
308 
309 
310
311                                                                                                                           
Note: See TracBrowser for help on using the repository browser.