chem/KPP/kpp/kpp-2.1/int/oldies/rodas3.c

   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
  12 void 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
  23 void 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
  53 int 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
 120 while(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