chem/KPP/kpp/kpp-2.1/int/oldies/ros3.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         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
  54 int ROS3(int N,KPP_REAL T,KPP_REAL Tnext,KPP_REAL Hmin,KPP_REAL Hmax,
  55 KPP_REAL Hstart,KPP_REAL y[NVAR],KPP_REAL AbsTol[NVAR],KPP_REAL RelTol[NVAR],
  56 int INFO[5],void (*forfun)(int,KPP_REAL,KPP_REAL [],KPP_REAL []) ,
  57 void (*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
 141 while( 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 /*
 289 this 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