source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/oldies/ros3.c @ 3423

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

Merge of branch palm4u into trunk

File size: 8.9 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        void (*forfun)(int,KPP_REAL,KPP_REAL [],KPP_REAL []);
10        void (*forjac)(int,KPP_REAL,KPP_REAL [],KPP_REAL []);
11
12
13
14    void FUNC_CHEM(int N,KPP_REAL T,KPP_REAL Y[NVAR],KPP_REAL P[NVAR])
15       {
16        KPP_REAL Told;
17        Told = TIME;
18        TIME = T;
19        Update_SUN();
20        Update_PHOTO();
21        Fun( Y, FIX, RCONST, P );
22        TIME = Told;
23       }/*  function fun ends here  */
24       
25
26    void JAC_CHEM(int N,KPP_REAL T,KPP_REAL Y[NVAR],KPP_REAL J[LU_NONZERO])
27       { 
28        KPP_REAL Told;
29        Told = TIME;
30        TIME = T;
31        Update_SUN();
32        Update_PHOTO();
33        Jac_SP( Y, FIX, RCONST, J );
34        TIME = Told;
35       } 
36       
37
38      INTEGRATE( KPP_REAL TIN, KPP_REAL TOUT )
39       {
40       
41        /*    TIN - Start Time          */
42        /*    TOUT - End Time           */
43   
44        int INFO[5];
45        forfun = &FUNC_CHEM;
46        forjac = &JAC_CHEM;     
47        INFO[0] = Autonomous;
48        ROS3(NVAR,TIN,TOUT,STEPMIN,STEPMAX,STEPMIN,VAR,ATOL,RTOL,INFO
49        ,forfun,forjac);
50       }        /* function integrate ends here  */
51     
52
53     
54int ROS3(int N,KPP_REAL T,KPP_REAL Tnext,KPP_REAL Hmin,KPP_REAL Hmax,
55KPP_REAL Hstart,KPP_REAL y[NVAR],KPP_REAL AbsTol[NVAR],KPP_REAL RelTol[NVAR],
56int INFO[5],void (*forfun)(int,KPP_REAL,KPP_REAL [],KPP_REAL []) ,
57void (*forjac)(int,KPP_REAL,KPP_REAL [],KPP_REAL []) )
58  {
59   
60/* 
61
62       L-stable Rosenbrock 3(2), with
63     strongly A-stable embedded formula for error control. 
64
65     All the arguments aggree with the KPP syntax.
66
67  INPUT ARGUMENTS:
68     y = Vector of (NVAR) concentrations, contains the
69         initial values on input
70     [T, Tnext] = the integration interval
71     Hmin, Hmax = lower and upper bounds for the selected step-size.
72          Note that for Step = Hmin the current computed
73          solution is unconditionally accepted by the error
74          control mechanism.
75     AbsTol, RelTol = (NVAR) dimensional vectors of
76          componentwise absolute and relative tolerances.
77     FUNC_CHEM = name of routine of derivatives. KPP syntax.
78          See the header below.
79     JAC_CHEM = name of routine that computes the Jacobian, in
80          sparse format. KPP syntax. See the header below.
81     Info(1) = 1  for  autonomous   system
82             = 0  for nonautonomous system
83
84  OUTPUT ARGUMENTS:
85     y = the values of concentrations at Tend.
86     T = equals Tend on output.
87     Info(2) = # of FUNC_CHEM calls.
88     Info(3) = # of JAC_CHEM calls.
89     Info(4) = # of accepted steps.
90     Info(5) = # of rejected steps.
91   
92     Adrian Sandu, April 1996
93     The Center for Global and Regional Environmental Research
94*/
95
96      KPP_REAL K1[NVAR], K2[NVAR], K3[NVAR];
97      KPP_REAL F1[NVAR], JAC[LU_NONZERO];
98      KPP_REAL ghinv,uround,dround,c43,x1,x2,x3,ytol;
99      KPP_REAL gam,c21,c31,c32,b1,b2,b3,d1,d2,d3,a21,a31,a32,alpha2,alpha3,
100                   g1,g2,g3;           
101      KPP_REAL ynew[NVAR];
102      KPP_REAL H, Hold, Tplus,tau;
103      KPP_REAL ERR, factor, facmax;
104      int n,nfcn,njac,Naccept,Nreject,i,j,ier;
105      char IsReject,Autonomous;
106
107
108/*     Initialization of counters, etc.         */
109      Autonomous = (INFO[0] == 1);         
110      uround = (double)1.e-15;
111      dround = DSQRT(uround);
112      H = MAX( (double)1.e-8, Hstart);
113      Tplus = T;
114      IsReject = 0;
115      Naccept  = 0;
116      Nreject  = 0;
117      nfcn     = 0;
118      njac     = 0;
119      gam  =   (double)  (.43586652150845899941601945119356e+00);
120      c21  =   (double) -(.10156171083877702091975600115545e+01);
121      c31  =   (double)  (.40759956452537699824805835358067e+01);
122      c32  =   (double)  (.92076794298330791242156818474003e+01);
123       b1  =   (double)  (.10000000000000000000000000000000e+01);
124       b2  =   (double)  (.61697947043828245592553615689730e+01);
125       b3  =   (double) -(.42772256543218573326238373806514e+00);
126       d1  =   (double)  (.50000000000000000000000000000000e+00);
127       d2  =   (double) -(.29079558716805469821718236208017e+01);
128       d3  =   (double)  (.22354069897811569627360909276199e+00);
129       a21 =   (double) 1.e0;
130       a31 =   (double) 1.e0;
131       a32 =   (double) 0.e0;
132       alpha2 = gam;
133       alpha3 = gam;
134       g1  =   (double)  (.43586652150845899941601945119356e+00);
135       g2  =   (double)  (.24291996454816804366592249683314e+00);
136       g3  =   (double)  (.21851380027664058511513169485832e+01);
137
138
139/* === Starting the time loop ===      */
140
141while( T < Tnext )
142 {
143  ten : 
144       Tplus = T + H;
145       
146       if ( Tplus > Tnext ) 
147       {
148          H = Tnext - T;
149          Tplus = Tnext;
150       }
151
152       (*forjac)(NVAR, T, y, JAC);
153       
154       njac = njac+1;
155       ghinv = (double)-1.0e0/(gam*H);
156
157       for(j=0;j<LU_NONZERO;j++)
158           JAC[j] = -JAC[j];
159       
160     
161
162       for(j=0;j<NVAR;j++)
163           JAC[LU_DIAG[j]] = JAC[LU_DIAG[j]] - ghinv;
164               
165       
166       
167       ier = KppDecomp (JAC);
168 
169       if ( ier != 0)
170       {
171        if( H > Hmin ) 
172        {
173         H = (double)5.0e-1*H; 
174         goto ten;         
175        }
176        else
177         {
178          printf("IER <> 0 , H = %d", H);
179         } 
180       }/* main ier if ends*/   
181               
182
183       (*forfun)(NVAR, T, y, F1);
184       
185
186/* ====== NONAUTONOMOUS CASE =============== */
187       if( Autonomous == 0 ) 
188       {
189         tau =(double) (dround*MAX( (double)1.0e-6, dabs(T) ) * signum(T) );
190
191         (*forfun)(NVAR, T+tau, y, K2);
192         
193         nfcn=nfcn+1;
194         
195         for(j=0;j<NVAR;j++)
196             K3[j] = ( K2[j]-F1[j] )/tau;
197           
198 
199
200/* ----- STAGE 1 (NONAUTONOMOUS) ----- */
201         x1 = (double)g1*H;
202         
203         for(j=0;j<NVAR;j++)
204             K1[j] =  F1[j] + x1*K3[j];
205           
206         KppSolve (JAC, K1);
207     
208/* ----- STAGE 2 (NONAUTONOMOUS) -----    */
209        for(j = 0;j<NVAR;j++)
210            ynew[j] = y[j] + K1[j]; 
211       
212       (*forfun)(NVAR, T+gam*H, ynew, F1);
213       nfcn=nfcn+1;
214       x1 = (double)(c21/H);
215       x2 = (double)(g2*H);
216       for(j = 0;j<NVAR;j++)
217           K2[j] = F1[j] + x1*K1[j] + x2*K3[j];
218       
219       KppSolve (JAC, K2);
220       
221/* ----- STAGE 3  (NONAUTONOMOUS) -----        */
222       x1 = (double)(c31/H);
223       x2 = (double)(c32/H);
224       x3 = (double)(g3*H);
225       for(j=0;j<NVAR;j++)
226           K3[j] = F1[j] + x1*K1[j] + x2*K2[j] + x3*K3[j];
227       
228       KppSolve (JAC, K3);
229      }/* "if" nonautonomous case ends here */ 
230
231
232
233/* ====== AUTONOMOUS CASE ===============  */
234
235       else
236       {
237       
238/* ----- STAGE 1 (AUTONOMOUS) -----   */
239
240         for(j = 0;j < NVAR;j++)
241             K1[j] =  F1[j]; 
242         
243         KppSolve (JAC, K1);
244 
245     
246/* ----- STAGE 2 (AUTONOMOUS) ----- */
247       for(j = 0;j < NVAR;j++)
248           ynew[j] = y[j] + K1[j];
249       
250       
251       (*forfun)(NVAR, T + gam*H, ynew, F1);
252       nfcn=nfcn+1;
253         x1 = (double)c21/H;
254         for(j = 0;j < NVAR;j++)
255             K2[j] = F1[j] + x1*K1[j]; 
256         
257         KppSolve (JAC, K2);
258       
259/* ----- STAGE 3  (AUTONOMOUS) -----   */
260
261       x1 = (double)(c31/H);
262       x2 = (double)(c32/H);
263       for(j = 0;j < NVAR;j++)
264           K3[j] = F1[j] + x1*K1[j] + x2*K2[j];
265       
266       KppSolve (JAC, K3);
267
268       }/*  Autonomousous case ends here  */
269       
270       
271    /*   ---- The Solution ---      */
272
273        for(j = 0;j < NVAR;j++)
274            ynew[j] = y[j] + b1*K1[j] + b2*K2[j] + b3*K3[j]; 
275       
276
277
278/* ====== Error estimation ========   */
279
280        ERR=(double)0.e0;
281        for(i=0;i<NVAR;i++)
282        {
283           ytol = AbsTol[i] + RelTol[i]*dabs(ynew[i]);
284           ERR = (double)(ERR+ pow( (double) ( (d1*K1[i]+d2*K2[i]+d3*K3[i])/ytol ) , 2 ));
285        }     
286        ERR = (double)MAX( uround, DSQRT( ERR/NVAR ) );
287
288/*
289this is the library i am linkin it to
290[sdmehra@herbert small_strato]$ ldd small_strato
291        libm.so.6 => /lib/libm.so.6 (0x40015000)
292        libc.so.6 => /lib/libc.so.6 (0x40032000)
293       /lib/ld-linux.so.2 => /lib/ld-linux.so.2 (0x40000000)                                                             
294*/
295
296/* ======= Choose the stepsize ===============================   */
297
298        factor = 0.9/pow( ERR , (1.e0/3.e0) );
299        if(IsReject == 1)
300           facmax = (double)1.0;
301           
302        else
303           facmax = (double)10.0;
304         
305           
306        factor = (double)MAX( 1.0e-1, MIN(factor,facmax) );
307        Hold = H;
308        H = (double)MIN( Hmax, MAX(Hmin,factor*H) );
309           
310
311/* ======= Rejected/Accepted Step ============================   */
312
313        if ( (ERR > 1) && (Hold > Hmin) ) 
314        {
315          IsReject = 1;
316          Nreject  = Nreject+1;
317        }
318        else
319        {
320          IsReject = 0;
321         
322         for(i = 0;i < NVAR;i++)
323             y[i]  = ynew[i];
324                     
325          T = Tplus;     
326          Naccept = Naccept+1;
327
328        }/*   else should end here */
329       
330
331/* ======= End of the time loop ===============================  */
332 
333   } /*      while loop (T < Tnext) ends here    */
334     
335     
336/* ======= Output Information =================================   */
337      INFO[1] = nfcn;
338      INFO[2] = njac;
339      INFO[3] = Naccept;
340      INFO[4] = Nreject;
341     
342    } /* function rodas ends here   */     
343
344 
Note: See TracBrowser for help on using the repository browser.