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