chem/KPP/kpp/kpp-2.1/int/oldies/ros2.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         #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
  17 void 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
  28 void 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
  60 int 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