[2696] | 1 | #define MAX(a,b) ( ((a) >= (b)) ?(a):(b) ) |
---|
| 2 | #define MIN(b,c) ( ((b) < (c)) ?(b):(c) ) |
---|
| 3 | #define ABS(x) ( ((x) >= 0 ) ?(x):(-x) ) |
---|
| 4 | #define SQRT(d) ( pow((d),0.5) ) |
---|
| 5 | /* SIGN transfer function */ |
---|
| 6 | #define SIGN(x,y) (((y) >= 0 ) ?(ABS(x)):(-ABS(x)) ) |
---|
| 7 | |
---|
| 8 | /*~~> Numerical constants */ |
---|
| 9 | #define ZERO (KPP_REAL)0.0 |
---|
| 10 | #define ONE (KPP_REAL)1.0 |
---|
| 11 | |
---|
| 12 | /*~~~> Statistics on the work performed by the Runge-Kutta method */ |
---|
| 13 | #define Nfun 0 |
---|
| 14 | #define Njac 1 |
---|
| 15 | #define Nstp 2 |
---|
| 16 | #define Nacc 3 |
---|
| 17 | #define Nrej 4 |
---|
| 18 | #define Ndec 5 |
---|
| 19 | #define Nsol 6 |
---|
| 20 | #define Nsng 7 |
---|
| 21 | #define Ntexit 0 |
---|
| 22 | #define Nhacc 1 |
---|
| 23 | #define Nhnew 2 |
---|
| 24 | |
---|
| 25 | /*~~~> Runge-Kutta method parameters */ |
---|
| 26 | #define RKmax 3 |
---|
| 27 | |
---|
| 28 | #define R2A 1 |
---|
| 29 | #define R1A 2 |
---|
| 30 | #define L3C 3 |
---|
| 31 | #define GAU 4 |
---|
| 32 | #define L3A 5 |
---|
| 33 | |
---|
| 34 | int rkMethod, |
---|
| 35 | SdirkError; |
---|
| 36 | KPP_REAL rkT[RKmax][RKmax], |
---|
| 37 | rkTinv[RKmax][RKmax], |
---|
| 38 | rkTinvAinv[RKmax][RKmax], |
---|
| 39 | rkAinvT[RKmax][RKmax], |
---|
| 40 | rkA[RKmax+1][RKmax+1], |
---|
| 41 | rkB[RKmax+1], |
---|
| 42 | rkC[RKmax+1], |
---|
| 43 | rkD[RKmax+1], |
---|
| 44 | rkE[RKmax+1], |
---|
| 45 | rkBgam[RKmax+2], |
---|
| 46 | rkBhat[RKmax+2], |
---|
| 47 | rkTheta[RKmax+1], |
---|
| 48 | rkF[RKmax+2], |
---|
| 49 | rkGamma, |
---|
| 50 | rkAlpha, |
---|
| 51 | rkBeta, |
---|
| 52 | rkELO; |
---|
| 53 | /*~~~> Function headers */ |
---|
| 54 | // void INTEGRATE(KPP_REAL TIN, KPP_REAL TOUT, int ICNTRL_U[], KPP_REAL RCNTRL_U[], |
---|
| 55 | // int ISTATUS_U[], KPP_REAL RSTATUS_U[], int IERR_U); |
---|
| 56 | void INTEGRATE(KPP_REAL TIN, KPP_REAL TOUT); |
---|
| 57 | void RungeKutta(int N, KPP_REAL T, KPP_REAL Tend, KPP_REAL Y[], |
---|
| 58 | KPP_REAL RelTol[], KPP_REAL AbsTol[], KPP_REAL RCNTRL[], |
---|
| 59 | int ICNTRL[], KPP_REAL RSTATUS[], int ISTATUS[], int* IERR); |
---|
| 60 | void RK_Integrator(int N, KPP_REAL* T, KPP_REAL Tend, KPP_REAL Y[], |
---|
| 61 | KPP_REAL AbsTol[], KPP_REAL RelTol[], int ITOL, |
---|
| 62 | int ISTATUS[], KPP_REAL RSTATUS[], KPP_REAL Hmin, |
---|
| 63 | KPP_REAL Hmax, KPP_REAL Hstart, KPP_REAL Roundoff, |
---|
| 64 | int Max_no_steps, int NewtonMaxit, int StartNewton, |
---|
| 65 | int Gustafsson, KPP_REAL ThetaMin, KPP_REAL NewtonTol, |
---|
| 66 | KPP_REAL FacSafe, KPP_REAL FacMax, KPP_REAL FacMin, |
---|
| 67 | KPP_REAL FacRej, KPP_REAL Qmin, KPP_REAL Qmax, int* IERR); |
---|
| 68 | void RK_ErrorMsg(int Code, KPP_REAL T, KPP_REAL H, int* IERR); |
---|
| 69 | void RK_ErrorScale(int N, int ITOL, KPP_REAL AbsTol[], KPP_REAL RelTol[], |
---|
| 70 | KPP_REAL Y[], KPP_REAL SCAL[]); |
---|
| 71 | /*void RK_Transform(int N, KPP_REAL Tr[][RKmax], KPP_REAL Z1[], KPP_REAL Z2[], |
---|
| 72 | KPP_REAL Z3[], KPP_REAL W1[], KPP_REAL W2[], KPP_REAL W3[]);*/ |
---|
| 73 | void RK_Interpolate(char action[], int N, KPP_REAL H, KPP_REAL Hold, |
---|
| 74 | KPP_REAL Z1[], KPP_REAL Z2[], KPP_REAL Z3[], KPP_REAL CONT[][RKmax]); |
---|
| 75 | void RK_PrepareRHS(int N, KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL FO[], |
---|
| 76 | KPP_REAL Z1[], KPP_REAL Z2[], KPP_REAL Z3[], KPP_REAL R1[], KPP_REAL R2[], |
---|
| 77 | KPP_REAL R3[]); |
---|
| 78 | void RK_Decomp(int N, KPP_REAL H, KPP_REAL FJAC[], KPP_REAL E1[], |
---|
| 79 | int IP1[], KPP_REAL E2R[], KPP_REAL E2I[], |
---|
| 80 | int IP2[], int* ISING, int ISTATUS[]); |
---|
| 81 | void RK_Solve(int N, KPP_REAL H, KPP_REAL E1[], int IP1[], KPP_REAL E2R[], |
---|
| 82 | KPP_REAL E2I[], int IP2[], KPP_REAL R1[], KPP_REAL R2[], KPP_REAL R3[], |
---|
| 83 | int ISTATUS[]); |
---|
| 84 | void RK_ErrorEstimate(int N, KPP_REAL H, KPP_REAL T, KPP_REAL Y[], |
---|
| 85 | KPP_REAL FO[], KPP_REAL E1[], int IP1[], KPP_REAL Z1[], |
---|
| 86 | KPP_REAL Z2[], KPP_REAL Z3[], KPP_REAL SCAL[], KPP_REAL* Err, |
---|
| 87 | int FirstStep, int Reject, int ISTATUS[]); |
---|
| 88 | void Radau2A_Coefficients(); |
---|
| 89 | void Lobatto3C_Coefficients (); |
---|
| 90 | void Gauss_Coefficients(); |
---|
| 91 | void Radau1A_Coefficients(); |
---|
| 92 | void Lobatto3A_Coefficients(); |
---|
| 93 | void FUN_CHEM(KPP_REAL T, KPP_REAL V[], KPP_REAL FCT[]); |
---|
| 94 | void JAC_CHEM(KPP_REAL T, KPP_REAL V[], KPP_REAL JF[]); |
---|
| 95 | KPP_REAL RK_ErrorNorm(int N, KPP_REAL SCAL[], KPP_REAL DY[]); |
---|
| 96 | void Fun(KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[]); |
---|
| 97 | void Jac_SP(KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[]); |
---|
| 98 | void WCOPY(int N, KPP_REAL X[], int incX, KPP_REAL Y[], int incY); |
---|
| 99 | void WADD(int N, KPP_REAL Y[], KPP_REAL Z[], KPP_REAL TMP[]); |
---|
| 100 | void WAXPY(int N, KPP_REAL Alpha, KPP_REAL X[], int incX, KPP_REAL Y[], int incY ); |
---|
| 101 | KPP_REAL WLAMCH( char C ); |
---|
| 102 | void KppSolveCmplxR(KPP_REAL JVSR[], KPP_REAL JVSI[], KPP_REAL XR[], KPP_REAL XI[]); |
---|
| 103 | int KppDecompCmplxR(KPP_REAL *JVSR, KPP_REAL *JVSI); |
---|
| 104 | void Set2Zero(int N, KPP_REAL A[]); |
---|
| 105 | int KppDecomp( KPP_REAL A[] ); |
---|
| 106 | void KppSolve ( KPP_REAL A[], KPP_REAL b[] ); |
---|
| 107 | void Update_SUN(); |
---|
| 108 | void Update_RCONST(); |
---|
| 109 | void Update_PHOTO(); |
---|
| 110 | |
---|
| 111 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 112 | //void INTEGRATE(KPP_REAL TIN, KPP_REAL TOUT, int ICNTRL_U[], KPP_REAL RCNTRL_U[], |
---|
| 113 | // int ISTATUS_U[], KPP_REAL RSTATUS_U[], int IERR_U ) |
---|
| 114 | void INTEGRATE(KPP_REAL TIN, KPP_REAL TOUT ) |
---|
| 115 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 116 | /* RungeKutta - Fully Implicit 3-stage Runge-Kutta methods based on: |
---|
| 117 | * Radau-2A quadrature (order 5) |
---|
| 118 | * Radau-1A quadrature (order 5) |
---|
| 119 | * Lobatto-3C quadrature (order 4) |
---|
| 120 | * Gauss quadrature (order 6) |
---|
| 121 | By default the code employs the KPP sparse linear algebra routines |
---|
| 122 | Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) |
---|
| 123 | |
---|
| 124 | (C) Adrian Sandu, August 2005 |
---|
| 125 | Virginia Polytechnic Institute and State University |
---|
| 126 | Contact: sandu@cs.vt.edu |
---|
| 127 | Revised by Philipp Miehe and Adrian Sandu, May 2006 |
---|
| 128 | F90 to C translation by Tinting Jiang and Don Jacob, July 2006 |
---|
| 129 | This implementation is part of KPP - the Kinetic PreProcessor |
---|
| 130 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 131 | { |
---|
| 132 | int IERR; |
---|
| 133 | KPP_REAL RCNTRL[20], |
---|
| 134 | RSTATUS[20], |
---|
| 135 | T1, |
---|
| 136 | T2; |
---|
| 137 | int ICNTRL[20], |
---|
| 138 | ISTATUS[20]; |
---|
| 139 | static int Ntotal = 0; /* for printing the number of steps */ |
---|
| 140 | |
---|
| 141 | int i; |
---|
| 142 | for ( i = 0; i < 20; i++ ) { |
---|
| 143 | RCNTRL[i] = ZERO; |
---|
| 144 | ICNTRL[i] = 0; |
---|
| 145 | } /* for */ |
---|
| 146 | |
---|
| 147 | /*~~~> fine-tune the integrator: */ |
---|
| 148 | ICNTRL[1] = 0; /* 0=vector tolerances, 1=scalar tolerances */ |
---|
| 149 | ICNTRL[4] = 8; /* Max no. of Newton iterations */ |
---|
| 150 | ICNTRL[5] = 0; /* Starting values for Newton are interpolated(0) or zero(1) */ |
---|
| 151 | ICNTRL[9] = 1; /* 0 - classic or 1 - SDIRK error estimation */ |
---|
| 152 | ICNTRL[10] = 0; /* Gustaffson(0) or classic(1) controller */ |
---|
| 153 | |
---|
| 154 | // /*~~~> if optional parameters are given, and if they are >0, |
---|
| 155 | // then use them to overwrite default settings */ |
---|
| 156 | // if (ICNTRL_U != NULL) { |
---|
| 157 | // for ( i = 0; i < 20; i++) { |
---|
| 158 | // if (ICNTRL_U[i] > 0) { |
---|
| 159 | // ICNTRL[i] = ICNTRL_U[i]; |
---|
| 160 | // } /* end if */ |
---|
| 161 | // } /* end for */ |
---|
| 162 | // } /* end if */ |
---|
| 163 | // if (RCNTRL_U != NULL) { |
---|
| 164 | // for (i = 0; i < 20; i++) { |
---|
| 165 | // if (RCNTRL_U[i] > 0) { |
---|
| 166 | // RCNTRL[i] = RCNTRL_U[i]; |
---|
| 167 | // } /* end if */ |
---|
| 168 | // } /* end for */ |
---|
| 169 | // } /* end if */ |
---|
| 170 | |
---|
| 171 | T1 = TIN; |
---|
| 172 | /*printf("T1=%f\n", T1);*/ |
---|
| 173 | T2 = TOUT; |
---|
| 174 | /*printf("T2=%f\n", T2);*/ |
---|
| 175 | RungeKutta(NVAR, T1, T2, VAR, RTOL, ATOL, RCNTRL,ICNTRL,RSTATUS,ISTATUS, &IERR); |
---|
| 176 | |
---|
| 177 | Ntotal += ISTATUS[Nstp]; |
---|
| 178 | printf("NSTEPS=%d (%d) O3=%E ", ISTATUS[Nstp], Ntotal, VAR[ind_O3]); |
---|
| 179 | |
---|
| 180 | // /* if optional parameters are given for output |
---|
| 181 | // * use them to store information in them */ |
---|
| 182 | // if (ISTATUS_U != NULL) { |
---|
| 183 | // for (i = 0; i < 20; i++) { |
---|
| 184 | // ISTATUS_U[i] = ISTATUS[i]; |
---|
| 185 | // } /* end for */ |
---|
| 186 | // } /* end if */ |
---|
| 187 | // if (RSTATUS_U != NULL) { |
---|
| 188 | // for (i = 0; i < 20; i++) { |
---|
| 189 | // RSTATUS_U[i] = RSTATUS[i]; |
---|
| 190 | // } /* end for */ |
---|
| 191 | // } /* end if */ |
---|
| 192 | // /*if (IERR_U != NULL) */ |
---|
| 193 | // IERR_U = IERR; |
---|
| 194 | |
---|
| 195 | if (IERR < 0) { |
---|
| 196 | printf("Runge-Kutta: Unsuccessful exit at T=%f(IERR=%d)", TIN, IERR); |
---|
| 197 | } /* end if */ |
---|
| 198 | |
---|
| 199 | } /* end INTEGRATE */ |
---|
| 200 | |
---|
| 201 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 202 | |
---|
| 203 | void RungeKutta(int N, KPP_REAL T, KPP_REAL Tend, KPP_REAL Y[], |
---|
| 204 | KPP_REAL RelTol[], KPP_REAL AbsTol[], KPP_REAL RCNTRL[], |
---|
| 205 | int ICNTRL[], KPP_REAL RSTATUS[], int ISTATUS[], int* IERR) |
---|
| 206 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 207 | |
---|
| 208 | This implementation is based on the book and the code Radau5: |
---|
| 209 | |
---|
| 210 | E. HAIRER AND G. WANNER |
---|
| 211 | "SOLVING ORDINARY DIFFERENTIAL EQUATIONS II. |
---|
| 212 | STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS." |
---|
| 213 | SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14, |
---|
| 214 | SPRINGER-VERLAG (1991) |
---|
| 215 | |
---|
| 216 | UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES |
---|
| 217 | CH-1211 GENEVE 24, SWITZERLAND |
---|
| 218 | E-MAIL: HAIRER@DIVSUN.UNIGE.CH, WANNER@DIVSUN.UNIGE.CH |
---|
| 219 | |
---|
| 220 | Methods: |
---|
| 221 | * Radau-2A quadrature (order 5) |
---|
| 222 | * Radau-1A quadrature (order 5) |
---|
| 223 | * Lobatto-3C quadrature (order 4) |
---|
| 224 | * Gauss quadrature (order 6) |
---|
| 225 | |
---|
| 226 | (C) Adrian Sandu, August 2005 |
---|
| 227 | Virginia Polytechnic Institute and State University |
---|
| 228 | Contact: sandu@cs.vt.edu |
---|
| 229 | This implementation is part of KPP - the Kinetic PreProcessor |
---|
| 230 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 231 | |
---|
| 232 | *~~~> INPUT ARGUMENTS: |
---|
| 233 | ---------------- |
---|
| 234 | Note: For input parameters equal to zero the default values of the |
---|
| 235 | corresponding variables are used. |
---|
| 236 | |
---|
| 237 | N Dimension of the system |
---|
| 238 | T Initial time value |
---|
| 239 | |
---|
| 240 | Tend Final T value (Tend-T may be positive or negative) |
---|
| 241 | |
---|
| 242 | Y(N) Initial values for Y |
---|
| 243 | |
---|
| 244 | RelTol, AbsTol Relative and absolute error tolerances. |
---|
| 245 | for ICNTRL[1] = 0: AbsTol, RelTol are N-dimensional vectors |
---|
| 246 | = 1: AbsTol, RelTol are scalars |
---|
| 247 | *~~~> Integer input parameters: |
---|
| 248 | |
---|
| 249 | ICNTRL[0] = not used |
---|
| 250 | |
---|
| 251 | ICNTRL[1] = 0: AbsTol, RelTol are NVAR-dimensional vectors |
---|
| 252 | = 1: AbsTol, RelTol are scalars |
---|
| 253 | |
---|
| 254 | ICNTRL[2] = RK method selection |
---|
| 255 | = 1: Radau-2A (the default) |
---|
| 256 | = 2: Radau-1A |
---|
| 257 | = 3: Lobatto-3C |
---|
| 258 | = 4: Gauss |
---|
| 259 | = 5: Lobatto-3A (not yet implemented) |
---|
| 260 | |
---|
| 261 | ICNTRL[3] -> maximum number of integration steps |
---|
| 262 | For ICNTRL[3] = 0 the default value of 10000 is used |
---|
| 263 | |
---|
| 264 | ICNTRL[4] -> maximum number of Newton iterations |
---|
| 265 | For ICNTRL[4] = 0 the default value of 8 is used |
---|
| 266 | |
---|
| 267 | ICNTRL[5] -> starting values of Newton iterations: |
---|
| 268 | ICNTRL[5] = 0: starting values are obtained from |
---|
| 269 | the extrapolated collocation solution |
---|
| 270 | (the default) |
---|
| 271 | ICNTRL[5] = 1: starting values are zero |
---|
| 272 | |
---|
| 273 | ICNTRL[9] -> switch for error estimation strategy |
---|
| 274 | ICNTRL[9] = 0: one additional stage at c = 0, |
---|
| 275 | see Hairer (default) |
---|
| 276 | ICNTRL[9] = 1: two additional stages at c = 0, |
---|
| 277 | and SDIRK at c = 1, stiffly accurate |
---|
| 278 | |
---|
| 279 | ICNTRL[10] -> switch for step size strategy |
---|
| 280 | ICNTRL[10] = 0: mod. predictive controller (Gustafsson, default) |
---|
| 281 | ICNTRL[10] = 1: classical step size control |
---|
| 282 | the choice 1 seems to produce safer results; |
---|
| 283 | for simple problems, the choice 2 produces |
---|
| 284 | often slightly faster runs |
---|
| 285 | |
---|
| 286 | *~~~> Real input parameters: |
---|
| 287 | RCNTRL[0] -> Hmin, lower bound for the integration step size |
---|
| 288 | (highly recommended to keep Hmin = ZERO, the default) |
---|
| 289 | |
---|
| 290 | RCNTRL[1] -> Hmax, upper bound for the integration step size |
---|
| 291 | |
---|
| 292 | RCNTRL[2] -> Hstart, the starting step size |
---|
| 293 | |
---|
| 294 | RCNTRL[3] -> FacMin, lower bound on step decrease factor |
---|
| 295 | (default=0.2) |
---|
| 296 | |
---|
| 297 | RCNTRL[4] -> FacMax, upper bound on step increase factor |
---|
| 298 | (default=6) |
---|
| 299 | |
---|
| 300 | RCNTRL[5] -> FacRej, step decrease factor after multiple rejections |
---|
| 301 | (default=0.1) |
---|
| 302 | |
---|
| 303 | RCNTRL[6] -> FacSafe, by which the new step is slightly smaller |
---|
| 304 | than the predicted value (default=0.9) |
---|
| 305 | |
---|
| 306 | RCNTRL[7] -> ThetaMin. If Newton convergence rate smaller |
---|
| 307 | than ThetaMin the Jacobian is not recomputed; |
---|
| 308 | (default=0.001) |
---|
| 309 | |
---|
| 310 | RCNTRL[8] -> NewtonTol, stopping criterion for Newton's method |
---|
| 311 | (default=0.03) |
---|
| 312 | |
---|
| 313 | RCNTRL[9] -> Qmin |
---|
| 314 | |
---|
| 315 | RCNTRL[10] -> Qmax. If Qmin < Hnew/Hold < Qmax, then the |
---|
| 316 | step size is kept constant and the LU factorization |
---|
| 317 | reused (default Qmin=1, Qmax=1.2) |
---|
| 318 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 319 | |
---|
| 320 | *~~~> OUTPUT ARGUMENTS: |
---|
| 321 | ----------------- |
---|
| 322 | T -> T value for which the solution has been computed |
---|
| 323 | (after successful return T=Tend). |
---|
| 324 | |
---|
| 325 | Y(N) -> Numerical solution at T |
---|
| 326 | |
---|
| 327 | Note: each call to RungeKutta adds the current no. of fcn calls |
---|
| 328 | to previous value of ISTATUS[0], and similar to other params. |
---|
| 329 | Set ISTATUS[1:9] = 0 before call to avoid this accumulation. |
---|
| 330 | |
---|
| 331 | ISTATUS[0] -> No. of function calls |
---|
| 332 | ISTATUS[1] -> No. of Jacobian calls |
---|
| 333 | ISTATUS[2] -> No. of steps |
---|
| 334 | ISTATUS[3] -> No. of accepted steps |
---|
| 335 | ISTATUS[4] -> No. of rejected steps (except at very beginning) |
---|
| 336 | ISTATUS[5] -> No. of LU decompositions |
---|
| 337 | ISTATUS[6] -> No. of forward/backward substitutions |
---|
| 338 | ISTATUS[7] -> No. of singular matrix decompositions |
---|
| 339 | |
---|
| 340 | RSTATUS[0] -> Texit, the time corresponding to the |
---|
| 341 | computed Y upon return |
---|
| 342 | RSTATUS[1] -> Hexit, last accepted step before exit |
---|
| 343 | RSTATUS[2] -> Hnew, last predicted step (not yet taken) |
---|
| 344 | For multiple restarts, use Hnew as Hstart |
---|
| 345 | in the subsequent run |
---|
| 346 | |
---|
| 347 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 348 | |
---|
| 349 | *~~~> RETURN VALUE (int): |
---|
| 350 | |
---|
| 351 | IERR -> Reports on successfulness upon return: |
---|
| 352 | = 1 for success |
---|
| 353 | < 0 for error (value equals error code) |
---|
| 354 | = -1 : Improper value for maximal no of steps |
---|
| 355 | = -2 : Improper value for maximal no of Newton iterations |
---|
| 356 | = -3 : Hmin/Hmax/Hstart must be positive |
---|
| 357 | = -4 : Improper values for FacMin/FacMax/FacSafe/FacRej |
---|
| 358 | = -5 : Improper value for ThetaMin |
---|
| 359 | = -6 : Newton stopping tolerance too small |
---|
| 360 | = -7 : Improper values for Qmin, Qmax |
---|
| 361 | = -8 : Tolerances are too small |
---|
| 362 | = -9 : No of steps exceeds maximum bound |
---|
| 363 | = -10 : Step size too small |
---|
| 364 | = -11 : Matrix is repeatedly singular |
---|
| 365 | = -12 : Non-convergence of Newton iterations |
---|
| 366 | = -13 : Requested RK method not implemented |
---|
| 367 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ |
---|
| 368 | { |
---|
| 369 | /*~~~> Control arguments */ |
---|
| 370 | /*printf("Starting RungeKutta\n");*/ |
---|
| 371 | int Max_no_steps; |
---|
| 372 | KPP_REAL Hmin, |
---|
| 373 | Hmax, |
---|
| 374 | Hstart, |
---|
| 375 | Qmin, |
---|
| 376 | Qmax, |
---|
| 377 | Roundoff, |
---|
| 378 | ThetaMin, |
---|
| 379 | NewtonTol, |
---|
| 380 | FacSafe, |
---|
| 381 | FacMin, |
---|
| 382 | FacMax, |
---|
| 383 | FacRej; |
---|
| 384 | /*~~~> Local variables */ |
---|
| 385 | int NewtonMaxit, |
---|
| 386 | ITOL, |
---|
| 387 | i, |
---|
| 388 | StartNewton, |
---|
| 389 | Gustafsson; |
---|
| 390 | |
---|
| 391 | *IERR = 0; |
---|
| 392 | for (i = 0; i < 20; i++) { |
---|
| 393 | ISTATUS[i] = 0; |
---|
| 394 | RSTATUS[i] = ZERO; |
---|
| 395 | } /* end for */ |
---|
| 396 | |
---|
| 397 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ |
---|
| 398 | |
---|
| 399 | /*~~~> ICNTRL[0] - autonomous system - not used */ |
---|
| 400 | /*~~~> ITOL: 1 for vector and 0 for scalar AbsTol/RelTol */ |
---|
| 401 | if (ICNTRL[1] == 0) { |
---|
| 402 | /*printf("Entering if ICNTRL[1] == 0\n");*/ |
---|
| 403 | ITOL = 1; |
---|
| 404 | } |
---|
| 405 | else { |
---|
| 406 | ITOL = 0; |
---|
| 407 | } /*end if */ |
---|
| 408 | /*~~~> Error control selection */ |
---|
| 409 | if (ICNTRL[9] == 0) { |
---|
| 410 | SdirkError = 0; |
---|
| 411 | } |
---|
| 412 | else { |
---|
| 413 | /*printf("Entering if ICNTRL[9] == 1\n");*/ |
---|
| 414 | SdirkError = 1; |
---|
| 415 | } /* end if */ |
---|
| 416 | /*~~~> Method Selection */ |
---|
| 417 | /*printf("Starting case ICNTRL[2]\n");*/ |
---|
| 418 | switch (ICNTRL[2]) { |
---|
| 419 | case 0: |
---|
| 420 | case 1: Radau2A_Coefficients(); /*printf("case 0 or 1\n");*/ |
---|
| 421 | break; |
---|
| 422 | case 2: Lobatto3C_Coefficients(); printf("case 2\n"); |
---|
| 423 | break; |
---|
| 424 | case 3: Gauss_Coefficients(); printf("case 3\n"); |
---|
| 425 | break; |
---|
| 426 | case 4: Radau1A_Coefficients(); printf("case 4\n"); |
---|
| 427 | break; |
---|
| 428 | case 5: Lobatto3A_Coefficients(); printf("case 5\n"); |
---|
| 429 | break; |
---|
| 430 | default: |
---|
| 431 | printf("\n ICNTRL[2]=%d\n", ICNTRL[2]); |
---|
| 432 | RK_ErrorMsg(-13, T, ZERO, IERR); |
---|
| 433 | } /* end switch */ |
---|
| 434 | |
---|
| 435 | /*~~~> Max_no_steps: the maximal number of time steps */ |
---|
| 436 | if (ICNTRL[3] == 0) { |
---|
| 437 | Max_no_steps = 200000; |
---|
| 438 | /*printf("Max_no_steps = %d\n", Max_no_steps);*/ |
---|
| 439 | } |
---|
| 440 | else { |
---|
| 441 | Max_no_steps = ICNTRL[3]; |
---|
| 442 | if (Max_no_steps <= 0) |
---|
| 443 | { |
---|
| 444 | printf("\n ICNTRL[3]=%d\n", ICNTRL[3]); |
---|
| 445 | RK_ErrorMsg(-1, T, ZERO, IERR); |
---|
| 446 | } /* end if */ |
---|
| 447 | } /* end if */ |
---|
| 448 | |
---|
| 449 | /*~~~> NewtonMaxit: maximal number of Newton iterations */ |
---|
| 450 | if (ICNTRL[4] == 0) |
---|
| 451 | NewtonMaxit = 8; |
---|
| 452 | else |
---|
| 453 | { |
---|
| 454 | NewtonMaxit = ICNTRL[4]; |
---|
| 455 | if (NewtonMaxit <= 0) { |
---|
| 456 | printf("\n ICNTRL[4]=%d\n", ICNTRL[4]); |
---|
| 457 | RK_ErrorMsg(-2, T, ZERO, IERR); |
---|
| 458 | } |
---|
| 459 | /*printf("NewtonMaxit = %d\n", NewtonMaxit);*/ |
---|
| 460 | } /* end if */ |
---|
| 461 | |
---|
| 462 | /*~~~> StartNewton: Use extrapolation for starting values of Newton iterations */ |
---|
| 463 | if (ICNTRL[5] == 0) { |
---|
| 464 | StartNewton = 1; |
---|
| 465 | /*printf("StartNewton = %d\n", StartNewton);*/ |
---|
| 466 | } |
---|
| 467 | else { |
---|
| 468 | StartNewton = 0; |
---|
| 469 | } /* end if */ |
---|
| 470 | |
---|
| 471 | /*~~~> Gustafsson: step size controller */ |
---|
| 472 | if (ICNTRL[10] == 0) { |
---|
| 473 | Gustafsson = 1; |
---|
| 474 | /*printf("Gustafsson = %d\n", Gustafsson);*/ |
---|
| 475 | } |
---|
| 476 | else { |
---|
| 477 | Gustafsson = 0; |
---|
| 478 | } /* end if */ |
---|
| 479 | |
---|
| 480 | /*~~~> Roundoff: smallest number s.t. 1.0 + Roundoff > 1.0 */ |
---|
| 481 | Roundoff = WLAMCH('E'); |
---|
| 482 | /*printf("Roundoff=%f\n", Roundoff);*/ |
---|
| 483 | /*~~~> Hmin = minimal step size */ |
---|
| 484 | if (RCNTRL[0] == ZERO) { |
---|
| 485 | /*printf("Entering if RCNTRL[0]=%f == ZERO\n", RCNTRL[0]);*/ |
---|
| 486 | Hmin = ZERO; |
---|
| 487 | } |
---|
| 488 | else { |
---|
| 489 | Hmin = MIN(ABS(RCNTRL[0]), ABS(Tend-T)); |
---|
| 490 | } /* end if */ |
---|
| 491 | /*~~~> Hmax = maximal step size */ |
---|
| 492 | if (RCNTRL[1] == ZERO) { |
---|
| 493 | /*printf("Entering if RCNTRL[1]=%f == ZERO\n", RCNTRL[1]);*/ |
---|
| 494 | Hmax = ABS(Tend-T); |
---|
| 495 | } |
---|
| 496 | else { |
---|
| 497 | Hmax = MIN(ABS(RCNTRL[1]),ABS(Tend-T)); |
---|
| 498 | } /* end if */ |
---|
| 499 | /*~~~> Hstart = starting step size */ |
---|
| 500 | if (RCNTRL[2] == ZERO) { |
---|
| 501 | /*printf("Entering if RCNTRL[2]=%f == ZERO\n", RCNTRL[2]);*/ |
---|
| 502 | Hstart = ZERO; |
---|
| 503 | } |
---|
| 504 | else { |
---|
| 505 | Hstart = MIN(ABS(RCNTRL[2]),ABS(Tend-T)); |
---|
| 506 | } /* end if */ |
---|
| 507 | /*~~~> FacMin: lower bound on step decrease factor */ |
---|
| 508 | if (RCNTRL[3] == ZERO) { |
---|
| 509 | /*printf("Entering if RCNTRL[3]=%f == ZERO\n", RCNTRL[3]);*/ |
---|
| 510 | FacMin = (KPP_REAL)0.2; |
---|
| 511 | } |
---|
| 512 | else { |
---|
| 513 | FacMin = RCNTRL[3]; |
---|
| 514 | } /* end if */ |
---|
| 515 | /*~~~> FacMax: upper bound on step increase factor */ |
---|
| 516 | if (RCNTRL[4] == ZERO) { |
---|
| 517 | /*printf("Entering if RCNTRL[4]=%f == ZERO\n", RCNTRL[4]);*/ |
---|
| 518 | FacMax = (KPP_REAL)8.0; |
---|
| 519 | } |
---|
| 520 | else { |
---|
| 521 | FacMax = RCNTRL[4]; |
---|
| 522 | } /* end if */ |
---|
| 523 | /*~~~> FacRej: step decrease factor after 2 consecutive rejections */ |
---|
| 524 | if (RCNTRL[5] == ZERO) { |
---|
| 525 | /*printf("Entering if RCNTRL[5]=%f == ZERO\n", RCNTRL[5]);*/ |
---|
| 526 | FacRej = (KPP_REAL)0.1; |
---|
| 527 | } |
---|
| 528 | else { |
---|
| 529 | FacRej = RCNTRL[5]; |
---|
| 530 | } /* end if */ |
---|
| 531 | /*~~~> FacSafe: by which the new step is slightly smaller |
---|
| 532 | than the predicted value */ |
---|
| 533 | if (RCNTRL[6] == ZERO) { |
---|
| 534 | /*printf("Entering if RCNTRL[6]=%f == ZERO\n", RCNTRL[6]);*/ |
---|
| 535 | FacSafe = (KPP_REAL)0.9; |
---|
| 536 | } |
---|
| 537 | else { |
---|
| 538 | FacSafe = RCNTRL[6]; |
---|
| 539 | } /* end if */ |
---|
| 540 | if ((FacMax < ONE) || (FacMin > ONE) || (FacSafe <= 1.0e-03) || (FacSafe >= ONE)) { |
---|
| 541 | printf("\n RCNTRL[3]=%f, RCNTRL[4]=%f, RCNTRL[5]=%f, RCNTRL[6]=%f\n", |
---|
| 542 | RCNTRL[3], RCNTRL[4], RCNTRL[5], RCNTRL[6]); |
---|
| 543 | RK_ErrorMsg(-4, T, ZERO, IERR); |
---|
| 544 | } /* end if */ |
---|
| 545 | |
---|
| 546 | /*~~~> ThetaMin: decides whether the Jacobian should be recomputed */ |
---|
| 547 | if (RCNTRL[7] == ZERO) |
---|
| 548 | ThetaMin = (KPP_REAL)1.0e-03; |
---|
| 549 | else { |
---|
| 550 | ThetaMin = RCNTRL[7]; |
---|
| 551 | if (ThetaMin <= (KPP_REAL)0.0 || ThetaMin >= (KPP_REAL)1.0) { |
---|
| 552 | printf("\n RCNTRL[7]=%f\n", RCNTRL[7]); |
---|
| 553 | RK_ErrorMsg(-5, T, ZERO, IERR); |
---|
| 554 | } |
---|
| 555 | } /* end if */ |
---|
| 556 | /*~~~> NewtonTol: stopping criterion for Newton's method */ |
---|
| 557 | if (RCNTRL[8] == ZERO) |
---|
| 558 | NewtonTol = (KPP_REAL)3.0e-02; |
---|
| 559 | else { |
---|
| 560 | NewtonTol = RCNTRL[8]; |
---|
| 561 | if (NewtonTol <= Roundoff) { |
---|
| 562 | printf("\n RCNTRL[8]=%f\n", RCNTRL[8]); |
---|
| 563 | RK_ErrorMsg(-6, T, ZERO, IERR); |
---|
| 564 | } |
---|
| 565 | } /* end if */ |
---|
| 566 | /*~~~> Qmin AND Qmax: IF Qmin < Hnew/Hold < Qmax then step size = const. */ |
---|
| 567 | if (RCNTRL[9] == ZERO) { |
---|
| 568 | Qmin = ONE; |
---|
| 569 | } |
---|
| 570 | else { |
---|
| 571 | Qmin = RCNTRL[9]; |
---|
| 572 | } /* end if */ |
---|
| 573 | if (RCNTRL[10] == ZERO) { |
---|
| 574 | Qmax = (KPP_REAL)1.2; |
---|
| 575 | } |
---|
| 576 | else { |
---|
| 577 | Qmax = RCNTRL[10]; |
---|
| 578 | } /* end if */ |
---|
| 579 | if (Qmin > ONE || Qmax < ONE) { |
---|
| 580 | printf("\n RCNTRL[9]=%f\n", Qmin); |
---|
| 581 | printf("\n RCNTRL[10]=%f\n", Qmax); |
---|
| 582 | RK_ErrorMsg(-7, T, ZERO, IERR); |
---|
| 583 | } /* end if */ |
---|
| 584 | /*~~~> Check if tolerances are reasonable */ |
---|
| 585 | if (ITOL == 0) { |
---|
| 586 | if ( AbsTol[0] <= ZERO || RelTol[0] <= ( ((KPP_REAL)10.0)*Roundoff) ) { |
---|
| 587 | printf("\n AbsTol=%f\n", AbsTol[0]); |
---|
| 588 | printf("\n RelTol=%f\n", RelTol[0]); |
---|
| 589 | RK_ErrorMsg(-8, T, ZERO, IERR); |
---|
| 590 | } |
---|
| 591 | } |
---|
| 592 | else { |
---|
| 593 | for (i = 0; i < N; i++) { |
---|
| 594 | if ( (AbsTol[i] <= ZERO) || (RelTol[i] <= ((KPP_REAL)10.0)*Roundoff) ) { |
---|
| 595 | printf("\n AbsTol[%d] = %f\n", i, AbsTol[i]); |
---|
| 596 | printf("\n AbsTol[%d] = %f\n", i, RelTol[i]); |
---|
| 597 | RK_ErrorMsg(-8, T, ZERO, IERR); |
---|
| 598 | } /* end if */ |
---|
| 599 | } /* end for */ |
---|
| 600 | } /* end if */ |
---|
| 601 | |
---|
| 602 | /*~~~> Parameters are wrong */ |
---|
| 603 | if (*IERR < 0) |
---|
| 604 | return; |
---|
| 605 | |
---|
| 606 | /*~~~> Call the core method */ |
---|
| 607 | RK_Integrator(N, &T, Tend, Y, AbsTol, RelTol, ITOL, ISTATUS, RSTATUS, |
---|
| 608 | Hmin, Hmax, Hstart, Roundoff, Max_no_steps, NewtonMaxit, |
---|
| 609 | StartNewton, Gustafsson, ThetaMin, NewtonTol, |
---|
| 610 | FacSafe, FacMax, FacMin, FacRej, Qmin, Qmax, IERR); |
---|
| 611 | /*printf("*IERR = %d\n", *IERR);*/ |
---|
| 612 | /*printf("Ending RungeKutta\n");*/ |
---|
| 613 | } /* RungeKutta */ |
---|
| 614 | |
---|
| 615 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 616 | void RK_Integrator( int N, |
---|
| 617 | /*~~~> Input: integration interval */ |
---|
| 618 | KPP_REAL* T, KPP_REAL Tend, KPP_REAL Y[], |
---|
| 619 | KPP_REAL AbsTol[], KPP_REAL RelTol[], int ITOL, |
---|
| 620 | /*~~~> Input: the initial condition at T; output: the solution at Tend */ |
---|
| 621 | int ISTATUS[], KPP_REAL RSTATUS[], KPP_REAL Hmin, KPP_REAL Hmax, |
---|
| 622 | KPP_REAL Hstart, KPP_REAL Roundoff, int Max_no_steps, int NewtonMaxit, |
---|
| 623 | int StartNewton, int Gustafsson, KPP_REAL ThetaMin, |
---|
| 624 | KPP_REAL NewtonTol, KPP_REAL FacSafe, KPP_REAL FacMax, KPP_REAL FacMin, |
---|
| 625 | KPP_REAL FacRej, KPP_REAL Qmin, KPP_REAL Qmax, int* IERR) |
---|
| 626 | { |
---|
| 627 | /*printf("Starting RK_Integrator\n");*/ |
---|
| 628 | KPP_REAL FJAC[LU_NONZERO], |
---|
| 629 | E1[LU_NONZERO], |
---|
| 630 | E2R[LU_NONZERO], |
---|
| 631 | E2I[LU_NONZERO]; |
---|
| 632 | KPP_REAL Z1[NVAR], |
---|
| 633 | Z2[NVAR], |
---|
| 634 | Z3[NVAR], |
---|
| 635 | Z4[NVAR], |
---|
| 636 | SCAL[NVAR], |
---|
| 637 | DZ1[NVAR], |
---|
| 638 | DZ2[NVAR], |
---|
| 639 | DZ3[NVAR], |
---|
| 640 | DZ4[NVAR], |
---|
| 641 | G[NVAR], |
---|
| 642 | TMP[NVAR], |
---|
| 643 | FO[NVAR]; |
---|
| 644 | KPP_REAL CONT[NVAR][RKmax], |
---|
| 645 | Tdirection, |
---|
| 646 | H, |
---|
| 647 | Hacc, |
---|
| 648 | Hnew, |
---|
| 649 | Hold, |
---|
| 650 | Fac, |
---|
| 651 | FacGus, |
---|
| 652 | Theta, |
---|
| 653 | Err, |
---|
| 654 | ErrOld, |
---|
| 655 | NewtonRate, |
---|
| 656 | NewtonIncrement, |
---|
| 657 | Hratio, |
---|
| 658 | Qnewton, |
---|
| 659 | NewtonPredictedErr, |
---|
| 660 | NewtonIncrementOld, |
---|
| 661 | ThetaSD; |
---|
| 662 | int IP1[NVAR], |
---|
| 663 | IP2[NVAR], |
---|
| 664 | NewtonIter, |
---|
| 665 | Nconsecutive, |
---|
| 666 | i; |
---|
| 667 | int ISING; |
---|
| 668 | int Reject, |
---|
| 669 | FirstStep, |
---|
| 670 | SkipJac, |
---|
| 671 | NewtonDone, |
---|
| 672 | SkipLU; |
---|
| 673 | |
---|
| 674 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 675 | |
---|
| 676 | /*~~~> INITIAL setting */ |
---|
| 677 | /*printf("ONE=%f Tend-*T=%f\n", ONE, Tend-*T);*/ |
---|
| 678 | Tdirection = SIGN(ONE, Tend-*T); |
---|
| 679 | /*printf("Tdirection=%f\n", Tdirection);*/ |
---|
| 680 | /*printf("Hmin=%g Hstart=%g \n", Hmin, Hstart);*/ |
---|
| 681 | /*printf("ABS(Hmin)=%g ABS(Hstart)=%g \n", ABS(Hmin), ABS(Hstart));*/ |
---|
| 682 | H = MIN( MAX(ABS(Hmin), ABS(Hstart)), Hmax ); |
---|
| 683 | if (ABS(H) <= ((KPP_REAL)10.0*Roundoff)) { |
---|
| 684 | /*printf("Entering if ABS(H)=%g <= 10.0*Roundoff=%g\n", ABS(H),(KPP_REAL)10.0*Roundoff);*/ |
---|
| 685 | H = (KPP_REAL)(1.0e-06); |
---|
| 686 | } /* end if */ |
---|
| 687 | H = SIGN(H, Tdirection); |
---|
| 688 | Hold = H; |
---|
| 689 | /*printf("Hold=%f\n", Hold);*/ |
---|
| 690 | Reject = 0; |
---|
| 691 | FirstStep = 1; |
---|
| 692 | SkipJac = 0; |
---|
| 693 | SkipLU = 0; |
---|
| 694 | if ((*T+H*((KPP_REAL)1.0001)-Tend)*Tdirection >= ZERO) { |
---|
| 695 | H = Tend - *T; |
---|
| 696 | } /* end if */ |
---|
| 697 | Nconsecutive = 0; |
---|
| 698 | RK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL); |
---|
| 699 | /*for(i=0; i<NVAR; i++) { |
---|
| 700 | printf("AbsTol=%g RelTol=%g Y=%g SCAL=%g \n", AbsTol[i], RelTol[i], Y[i], SCAL[i] ); |
---|
| 701 | }*/ |
---|
| 702 | |
---|
| 703 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 704 | /*~~~> Time loop begins */ |
---|
| 705 | /* while Tloop */ |
---|
| 706 | Tloop: while ( (Tend-*T)*Tdirection - Roundoff > ZERO ) { |
---|
| 707 | /*printf("Starting Tloop: (Tend-*T)*Tdirection - Roundoff=%g > ZERO\n", (Tend-*T)*Tdirection-Roundoff);*/ |
---|
| 708 | /*if ( Reject == 0 ) { */ |
---|
| 709 | FUN_CHEM(*T,Y,FO); |
---|
| 710 | ISTATUS[Nfun]++; |
---|
| 711 | /* } * end if */ |
---|
| 712 | if ( SkipLU == 0 ) { /* This time around skip the Jac update and LU */ |
---|
| 713 | /*~~~> Compute the Jacobian matrix */ |
---|
| 714 | /*printf("Entering if SkipLU == 0\n");*/ |
---|
| 715 | if ( SkipJac == 0 ) { |
---|
| 716 | /*printf("Entering if SkipJac == 0\n");*/ |
---|
| 717 | JAC_CHEM(*T,Y,FJAC); |
---|
| 718 | ISTATUS[Njac]++; |
---|
| 719 | } /* end if */ |
---|
| 720 | /*~~~> Compute the matrices E1 and E2 and their decompositions */ |
---|
| 721 | RK_Decomp(N,H,FJAC,E1,IP1,E2R,E2I,IP2,&ISING,ISTATUS); |
---|
| 722 | /*printf("ISING=%d\n", ISING);*/ |
---|
| 723 | if ( ISING != 0 ) { |
---|
| 724 | /*printf("Entering if ISING != 0\n");*/ |
---|
| 725 | ISTATUS[Nsng]++; |
---|
| 726 | Nconsecutive++; |
---|
| 727 | if (Nconsecutive >= 5) { |
---|
| 728 | RK_ErrorMsg(-12,*T,H,IERR); |
---|
| 729 | } |
---|
| 730 | H = H * ((KPP_REAL)0.5); |
---|
| 731 | Reject = 1; |
---|
| 732 | SkipJac = 1; |
---|
| 733 | SkipLU = 0; |
---|
| 734 | goto Tloop; |
---|
| 735 | } |
---|
| 736 | else { |
---|
| 737 | /*printf("Entering if ISING == 0\n");*/ |
---|
| 738 | Nconsecutive = 0; |
---|
| 739 | } /* end if */ |
---|
| 740 | } /* end if !SkipLU */ |
---|
| 741 | |
---|
| 742 | /*printf("NSTEPS=%d\n", ISTATUS[Nstp]);*/ |
---|
| 743 | ISTATUS[Nstp]++; |
---|
| 744 | /*printf("NSTEPS=%d\n", ISTATUS[Nstp]);*/ |
---|
| 745 | if (ISTATUS[Nstp] > Max_no_steps) { |
---|
| 746 | printf("\n Max number of time steps is = %d", Max_no_steps); |
---|
| 747 | RK_ErrorMsg(-9,*T,H,IERR); |
---|
| 748 | } /* end if */ |
---|
| 749 | if (((KPP_REAL)0.1)*ABS(H) <= ABS(*T)*Roundoff) |
---|
| 750 | RK_ErrorMsg(-10,*T,H,IERR); |
---|
| 751 | |
---|
| 752 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 753 | /*~~~> Loop for the simplified Newton iterations */ |
---|
| 754 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 755 | |
---|
| 756 | /*~~~> Starting values for Newton iteration */ |
---|
| 757 | if ( (FirstStep == 1) || (StartNewton == 0) ) { |
---|
| 758 | /*printf("Entering if FirstStep or Not StartNewton\n");*/ |
---|
| 759 | Set2Zero(N,Z1); |
---|
| 760 | Set2Zero(N,Z2); |
---|
| 761 | Set2Zero(N,Z3); |
---|
| 762 | } |
---|
| 763 | else |
---|
| 764 | { |
---|
| 765 | /* Evaluate quadratic polynomial */ |
---|
| 766 | RK_Interpolate("eval",N,H,Hold,Z1,Z2,Z3,CONT); |
---|
| 767 | } /* end if */ |
---|
| 768 | |
---|
| 769 | /*~~~> Initializations for Newton iteration */ |
---|
| 770 | NewtonDone = 0; |
---|
| 771 | Fac = (KPP_REAL)0.5; /* Step reduction if too many iterations */ |
---|
| 772 | |
---|
| 773 | /* for NewtonLoop */ |
---|
| 774 | for (NewtonIter = 1; NewtonIter <= NewtonMaxit; NewtonIter++) |
---|
| 775 | { |
---|
| 776 | /*printf("Starting NewtonIter loop: NewtonIter=%d\n", NewtonIter);*/ |
---|
| 777 | /* Prepare the right-hand side */ |
---|
| 778 | RK_PrepareRHS(N,*T,H,Y,FO,Z1,Z2,Z3,DZ1,DZ2,DZ3); |
---|
| 779 | |
---|
| 780 | /* Solve the linear systems */ |
---|
| 781 | RK_Solve(N,H,E1,IP1,E2R,E2I,IP2,DZ1,DZ2,DZ3,ISTATUS); |
---|
| 782 | |
---|
| 783 | NewtonIncrement = |
---|
| 784 | SQRT((RK_ErrorNorm(N,SCAL,DZ1) * RK_ErrorNorm(N,SCAL,DZ1) |
---|
| 785 | + RK_ErrorNorm(N,SCAL,DZ2) * RK_ErrorNorm(N,SCAL,DZ2) |
---|
| 786 | + RK_ErrorNorm(N,SCAL,DZ3) * RK_ErrorNorm(N,SCAL,DZ3)) |
---|
| 787 | / (KPP_REAL)3.0 ); |
---|
| 788 | /*printf( "NewtonIncrement=%g \n", NewtonIncrement );*/ |
---|
| 789 | if (NewtonIter == 1) { |
---|
| 790 | /*printf("Entering NewtonIter=%d == 1\n", NewtonIter);*/ |
---|
| 791 | Theta = ABS(ThetaMin); |
---|
| 792 | NewtonRate = (KPP_REAL)2.0; |
---|
| 793 | } |
---|
| 794 | else { |
---|
| 795 | /*printf("Entering else from NewtonIter=%d == 1\n", NewtonIter);*/ |
---|
| 796 | Theta = NewtonIncrement / NewtonIncrementOld; |
---|
| 797 | /*printf("Theta=%f\n", Theta);*/ |
---|
| 798 | if (Theta < (KPP_REAL)0.99) { |
---|
| 799 | /*printf( "Entering if Theta=%g < 0.99\n", Theta );*/ |
---|
| 800 | NewtonRate = Theta / (ONE-Theta); |
---|
| 801 | } |
---|
| 802 | else { /* Non-convergence of Newton: Theta too large */ |
---|
| 803 | break; /* EXIT NewtonLoop */ |
---|
| 804 | } /* end if */ |
---|
| 805 | if (NewtonIter < NewtonMaxit) { |
---|
| 806 | /*printf("Entering if NewtonIter=%d < NewtonMaxit=%d\n", NewtonIter, NewtonMaxit);*/ |
---|
| 807 | /* Predict error at the end of Newton process */ |
---|
| 808 | NewtonPredictedErr = NewtonIncrement |
---|
| 809 | *pow(Theta,(NewtonMaxit-NewtonIter))/(ONE-Theta); |
---|
| 810 | /*printf("NewtonRate=%g NewtonPredictedErr=%g \n", NewtonRate, NewtonPredictedErr);*/ |
---|
| 811 | if (NewtonPredictedErr >= NewtonTol) { |
---|
| 812 | /*printf( "Entering if NewtonPredictedErr=%g >= NewtonTol%g \n", NewtonPredictedErr, NewtonTol );*/ |
---|
| 813 | /* Non-convergence of Newton: predicted error too large */ |
---|
| 814 | Qnewton = MIN((KPP_REAL)10.0, NewtonPredictedErr/NewtonTol); |
---|
| 815 | Fac = (KPP_REAL)0.8*pow(Qnewton,(-ONE/(1+NewtonMaxit-NewtonIter))); |
---|
| 816 | break; /* EXIT NewtonLoop */ |
---|
| 817 | } /* end if */ |
---|
| 818 | } /* end if */ |
---|
| 819 | } /* end if */ |
---|
| 820 | |
---|
| 821 | NewtonIncrementOld = MAX(NewtonIncrement, Roundoff); |
---|
| 822 | /*printf("NewtonIncrementOld=%f\n", NewtonIncrementOld);*/ |
---|
| 823 | /*~~~> Update solution */ |
---|
| 824 | WAXPY(N,-ONE,DZ1,1,Z1,1); /* Z1 <- Z1 - DZ1 */ |
---|
| 825 | WAXPY(N,-ONE,DZ2,1,Z2,1); /* Z2 <- Z2 - DZ2 */ |
---|
| 826 | WAXPY(N,-ONE,DZ3,1,Z3,1); /* Z3 <- Z3 - DZ3 */ |
---|
| 827 | /*for(i=0; i<N; i++) |
---|
| 828 | printf("DZ1[%d]=%g Z1[%d]=%g\n", i, DZ1[i], i, Z1[i]);*/ |
---|
| 829 | /*~~~> Check error in Newton iterations */ |
---|
| 830 | /*printf( "NewtonRate=%g NewtonIncrement=%g NewtonTol=%g \n", NewtonRate, NewtonIncrement, NewtonTol );*/ |
---|
| 831 | NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol); |
---|
| 832 | /*printf( "NewtonDone=%d \n", NewtonDone );*/ |
---|
| 833 | if (NewtonDone == 1) |
---|
| 834 | break; /* Exit NewtonLoop */ |
---|
| 835 | if (NewtonIter == NewtonMaxit) |
---|
| 836 | { |
---|
| 837 | printf("\n Slow or no convergence in Newton Iteration:"); |
---|
| 838 | printf(" Max no. of Newton iterations reached"); |
---|
| 839 | } /* end if */ |
---|
| 840 | }/* end for NewtonLoop */ |
---|
| 841 | |
---|
| 842 | if ( NewtonDone == 0) |
---|
| 843 | { |
---|
| 844 | /*printf( "Entering if NewtonDone == 0\n" );*/ |
---|
| 845 | /*RK_ErrorMsg(-12,*T,H,IERR); */ |
---|
| 846 | H = Fac*H; |
---|
| 847 | Reject = 1; |
---|
| 848 | SkipJac = 1; |
---|
| 849 | SkipLU = 0; |
---|
| 850 | goto Tloop; |
---|
| 851 | } /* end if */ |
---|
| 852 | |
---|
| 853 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 854 | /*~~~> SDIRK Stage */ |
---|
| 855 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 856 | if (SdirkError == 1) |
---|
| 857 | { |
---|
| 858 | /*printf("Entering if SdirkError=%d==1\n", SdirkError);*/ |
---|
| 859 | /*~~~> Starting values for Newton iterations */ |
---|
| 860 | for (i = 0; i < N; i++) |
---|
| 861 | { |
---|
| 862 | Z4[i] = Z3[i]; |
---|
| 863 | } |
---|
| 864 | /*~~~> Prepare the loop-independent part of the right-hand side */ |
---|
| 865 | /* G = H*rkBgam(0)*FO + rkTheta(1)*Z1 |
---|
| 866 | + rkTheta(2)*Z2 + rkTheta(3)*Z3; */ |
---|
| 867 | Set2Zero(N,G); |
---|
| 868 | if (rkMethod != L3A) |
---|
| 869 | WAXPY(N,rkBgam[0]*H, FO,1,G,1); |
---|
| 870 | WAXPY(N,rkTheta[0],Z1,1,G,1); |
---|
| 871 | WAXPY(N,rkTheta[1],Z2,1,G,1); |
---|
| 872 | WAXPY(N,rkTheta[2],Z3,1,G,1); |
---|
| 873 | |
---|
| 874 | /*~~~> Initializations for Newton iteration */ |
---|
| 875 | NewtonDone = 0; |
---|
| 876 | Fac = (KPP_REAL)0.5; /* Step reduction factor if too many iterations */ |
---|
| 877 | |
---|
| 878 | /* for SDNewtonLoop */ |
---|
| 879 | for (NewtonIter = 1; NewtonIter <= NewtonMaxit; NewtonIter++) |
---|
| 880 | { |
---|
| 881 | /*printf("Starting NewtonIter loop: NewtonIter=%d\n", NewtonIter);*/ |
---|
| 882 | /*~~~> Prepare the loop-dependent part of right-hand side */ |
---|
| 883 | WADD(N,Y,Z4,TMP); /* TMP <- Y + Z4 */ |
---|
| 884 | FUN_CHEM(*T+H,TMP,DZ4); /* DZ4 <- Fun(Y+Z4) */ |
---|
| 885 | ISTATUS[Nfun]++; |
---|
| 886 | /* DZ4[1,N] = (D[1, N]-Z4[1,N])*(rkGamma/H) + DZ4[1,N]; */ |
---|
| 887 | WAXPY(N,-ONE*rkGamma/H,Z4,1,DZ4,1); |
---|
| 888 | WAXPY(N,rkGamma/H,G,1,DZ4,1); |
---|
| 889 | |
---|
| 890 | /*~~~> Solve the linear system */ |
---|
| 891 | KppSolve(E1, DZ4); |
---|
| 892 | /*~~~> Note: for a full matrix use Lapack: |
---|
| 893 | DGETRS('N', 5, 1, E1, N, IP1, DZ4, 5, ISING) */ |
---|
| 894 | |
---|
| 895 | /*~~~> Check convergence of Newton iterations */ |
---|
| 896 | NewtonIncrement = RK_ErrorNorm(N,SCAL,DZ4); |
---|
| 897 | /*printf( "NewtonIncrement=%g \n", NewtonIncrement );*/ |
---|
| 898 | if (NewtonIter == 1) { |
---|
| 899 | /*printf("Entering NewtonIter=%d == 1\n", NewtonIter);*/ |
---|
| 900 | ThetaSD = ABS(ThetaMin); |
---|
| 901 | NewtonRate = (KPP_REAL)2.0; |
---|
| 902 | } |
---|
| 903 | else { |
---|
| 904 | /*printf("Entering else from NewtonIter=%d == 1\n", NewtonIter);*/ |
---|
| 905 | ThetaSD = NewtonIncrement / NewtonIncrementOld; |
---|
| 906 | if (ThetaSD < (KPP_REAL)0.99) { |
---|
| 907 | /*printf( "Entering if Theta=%g < 0.99\n", Theta );*/ |
---|
| 908 | NewtonRate = ThetaSD / (ONE-ThetaSD); |
---|
| 909 | /* Predict error at the end of Newton process */ |
---|
| 910 | NewtonPredictedErr = NewtonIncrement |
---|
| 911 | *pow(ThetaSD,(NewtonMaxit-NewtonIter))/(ONE-ThetaSD); |
---|
| 912 | /*printf("NewtonRate=%g NewtonPredictedErr=%g \n", NewtonRate, NewtonPredictedErr);*/ |
---|
| 913 | if (NewtonPredictedErr >= NewtonTol) { |
---|
| 914 | /*printf( "Entering if NewtonPredictedErr=%g >= NewtonTol%g \n", NewtonPredictedErr, NewtonTol );*/ |
---|
| 915 | /* Non-convergence of Newton: predicted error too large */ |
---|
| 916 | /* printf("\n Error too large", NewtonPredictedErr); */ |
---|
| 917 | Qnewton = MIN((KPP_REAL)10.0,NewtonPredictedErr/NewtonTol); |
---|
| 918 | Fac = (KPP_REAL)0.8*pow(Qnewton,(-ONE/(1+NewtonMaxit-NewtonIter))); |
---|
| 919 | break; /* EXIT SDNewtonLoop */ |
---|
| 920 | } /* end if */ |
---|
| 921 | } |
---|
| 922 | /* Non-convergence of Newton: Theta too large */ |
---|
| 923 | else { |
---|
| 924 | /* prinf("\n Theta too large", ThetaSD); */ |
---|
| 925 | break; /* EXIT SDNewtonLoop */ |
---|
| 926 | } /* end if */ |
---|
| 927 | } /* end if */ |
---|
| 928 | NewtonIncrementOld = NewtonIncrement; |
---|
| 929 | /*printf("NewtonIncrementOld=%f\n", NewtonIncrementOld);*/ |
---|
| 930 | /* Update solution: Z4 <-- Z4 + DZ4; */ |
---|
| 931 | WAXPY(N,ONE,DZ4,1,Z4,1); |
---|
| 932 | |
---|
| 933 | /* Check error in Newton iterations */ |
---|
| 934 | NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol); |
---|
| 935 | /*printf( "NewtonDone=%d \n", NewtonDone );*/ |
---|
| 936 | if (NewtonDone == 1) |
---|
| 937 | break; /* EXIT SDNewtonLoop */ |
---|
| 938 | } /* end for SDNewtonLoop */ |
---|
| 939 | |
---|
| 940 | if ( NewtonDone == 0 ) { |
---|
| 941 | /*printf( "Entering if NewtonDone == 0\n" );*/ |
---|
| 942 | H = Fac*H; |
---|
| 943 | Reject = 1; |
---|
| 944 | SkipJac = 1; |
---|
| 945 | SkipLU = 0; |
---|
| 946 | goto Tloop; |
---|
| 947 | } /* end if */ |
---|
| 948 | } /* end if */ |
---|
| 949 | /*~~~> End of implified SDIRK Newton iterations */ |
---|
| 950 | |
---|
| 951 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 952 | /*~~~> Error estimation */ |
---|
| 953 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 954 | if (SdirkError == 1) { |
---|
| 955 | Set2Zero(N, DZ4); |
---|
| 956 | if (rkMethod == L3A) { |
---|
| 957 | for (i = 0; i < N; i++) |
---|
| 958 | DZ4[i] = H*rkF[0]*FO[i]; |
---|
| 959 | if (rkF[0] != ZERO) |
---|
| 960 | WAXPY(N, rkF[0], Z1, 1, DZ4, 1); |
---|
| 961 | if (rkF[1] != ZERO) |
---|
| 962 | WAXPY(N, rkF[1], Z2, 1, DZ4, 1); |
---|
| 963 | if (rkF[2] != ZERO) |
---|
| 964 | WAXPY(N, rkF[2], Z3, 1, DZ4, 1); |
---|
| 965 | for (i = 0; i < N; i++) |
---|
| 966 | TMP[i] = Y[i] + Z4[i]; |
---|
| 967 | FUN_CHEM(*T+H, TMP, DZ1); |
---|
| 968 | WAXPY(N, H*rkBgam[4], DZ1, 1, DZ4, 1); |
---|
| 969 | } |
---|
| 970 | else { |
---|
| 971 | /* DZ4(1,N) = rkD(1)*Z1 + rkD(2)*Z2 + rkD(3)*Z3 - Z4; */ |
---|
| 972 | if (rkD[0] != ZERO) |
---|
| 973 | WAXPY(N, rkD[0], Z1, 1, DZ4, 1); |
---|
| 974 | if (rkD[1] != ZERO) |
---|
| 975 | WAXPY(N, rkD[1], Z2, 1, DZ4, 1); |
---|
| 976 | if (rkD[2] != ZERO) |
---|
| 977 | WAXPY(N, rkD[2], Z3, 1, DZ4, 1); |
---|
| 978 | WAXPY(N, -ONE, Z4, 1, DZ4, 1); |
---|
| 979 | } /* end if */ |
---|
| 980 | Err = RK_ErrorNorm(N,SCAL,DZ4); |
---|
| 981 | } |
---|
| 982 | else |
---|
| 983 | { |
---|
| 984 | RK_ErrorEstimate(N,H,*T,Y,FO,E1,IP1,Z1,Z2,Z3,SCAL,&Err, |
---|
| 985 | FirstStep,Reject,ISTATUS); |
---|
| 986 | } /* end if */ |
---|
| 987 | |
---|
| 988 | /*~~~> Computation of new step size Hnew */ |
---|
| 989 | Fac = pow(Err, (-ONE/rkELO)) |
---|
| 990 | *MIN(FacSafe,(ONE+2*NewtonMaxit)/(NewtonIter+2*NewtonMaxit)); |
---|
| 991 | Fac = MIN(FacMax,MAX(FacMin,Fac)); |
---|
| 992 | Hnew = Fac*H; |
---|
| 993 | |
---|
| 994 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 995 | /*~~~> Accept/reject step */ |
---|
| 996 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 997 | // accept: |
---|
| 998 | if (Err < ONE) { /*~~~> STEP IS ACCEPTED */ |
---|
| 999 | FirstStep = 0; |
---|
| 1000 | ISTATUS[Nacc]++; |
---|
| 1001 | if (Gustafsson == 1) { |
---|
| 1002 | /*~~~> Predictive controller of Gustafsson */ |
---|
| 1003 | if (ISTATUS[Nacc] > 1) { |
---|
| 1004 | FacGus = FacSafe*(H/Hacc)*pow(Err*Err/ErrOld,(KPP_REAL)(-0.25)); |
---|
| 1005 | FacGus = MIN(FacMax,MAX(FacMin,FacGus)); |
---|
| 1006 | Fac = MIN(Fac,FacGus); |
---|
| 1007 | Hnew = Fac*H; |
---|
| 1008 | } /* end if */ |
---|
| 1009 | Hacc = H; |
---|
| 1010 | ErrOld = MAX((KPP_REAL)1.0e-02,Err); |
---|
| 1011 | } /* end if */ |
---|
| 1012 | Hold = H; |
---|
| 1013 | *T = *T + H; |
---|
| 1014 | /* Update solution: Y <- Y + sum(d_i Z-i) */ |
---|
| 1015 | if (rkD[0] != ZERO) |
---|
| 1016 | WAXPY(N, rkD[0], Z1, 1, Y, 1); |
---|
| 1017 | if (rkD[1] != ZERO) |
---|
| 1018 | WAXPY(N, rkD[1], Z2, 1, Y, 1); |
---|
| 1019 | if (rkD[2] != ZERO) |
---|
| 1020 | WAXPY(N, rkD[2], Z3, 1, Y, 1); |
---|
| 1021 | /* Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3 */ |
---|
| 1022 | if (StartNewton == 1) |
---|
| 1023 | RK_Interpolate("make",N,H,Hold,Z1,Z2,Z3,CONT); |
---|
| 1024 | RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL); |
---|
| 1025 | RSTATUS[Ntexit] = *T; |
---|
| 1026 | RSTATUS[Nhnew] = Hnew; |
---|
| 1027 | RSTATUS[Nhacc] = H; |
---|
| 1028 | Hnew = Tdirection*MIN( MAX(ABS(Hnew),Hmin), Hmax ); |
---|
| 1029 | if (Reject == 1) |
---|
| 1030 | Hnew = Tdirection*MIN(ABS(Hnew),ABS(H)); |
---|
| 1031 | Reject = 0; |
---|
| 1032 | if ((*T+Hnew/Qmin-Tend)*Tdirection >= ZERO) |
---|
| 1033 | H = Tend - *T; |
---|
| 1034 | else { |
---|
| 1035 | Hratio = Hnew/H; |
---|
| 1036 | /* Reuse the LU decomposition */ |
---|
| 1037 | SkipLU = (Theta <= ThetaMin) && (Hratio >= Qmin) && (Hratio <= Qmax); |
---|
| 1038 | if ( SkipLU == 0) |
---|
| 1039 | H = Hnew; |
---|
| 1040 | } /* end if */ |
---|
| 1041 | /* If convergence is fast enough, do not update Jacobian */ |
---|
| 1042 | /* SkipJac = (Theta <= ThetaMin); */ |
---|
| 1043 | SkipJac = 0; |
---|
| 1044 | } |
---|
| 1045 | |
---|
| 1046 | /*~~~> Step is rejected */ |
---|
| 1047 | else { |
---|
| 1048 | if ((FirstStep == 1) || (Reject == 1)) |
---|
| 1049 | H = FacRej*H; |
---|
| 1050 | else |
---|
| 1051 | H = Hnew; |
---|
| 1052 | Reject = 1; |
---|
| 1053 | SkipJac = 1; /* Skip if rejected - Jac is independent of H */ |
---|
| 1054 | SkipLU = 0; |
---|
| 1055 | if (ISTATUS[Nacc] >= 1) |
---|
| 1056 | ISTATUS[Nrej]++; |
---|
| 1057 | } /* end if accept */ |
---|
| 1058 | } /* while: time Tloop */ |
---|
| 1059 | |
---|
| 1060 | /*~~~> Successful exit */ |
---|
| 1061 | *IERR = 1; |
---|
| 1062 | |
---|
| 1063 | } /* RK_Integrator */ |
---|
| 1064 | |
---|
| 1065 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1066 | Handles all error messages and returns IERR = error Code |
---|
| 1067 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1068 | void RK_ErrorMsg(int Code, KPP_REAL T, KPP_REAL H, int* IERR) |
---|
| 1069 | { |
---|
| 1070 | Code = *IERR; |
---|
| 1071 | printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"); |
---|
| 1072 | printf("\nForced exit from RungeKutta due to the following error:\n"); |
---|
| 1073 | |
---|
| 1074 | switch (Code) { |
---|
| 1075 | case -1: |
---|
| 1076 | printf("--> Improper value for maximal no of steps"); |
---|
| 1077 | break; |
---|
| 1078 | case -2: |
---|
| 1079 | printf("--> Improper value for maximal no of Newton iterations"); |
---|
| 1080 | break; |
---|
| 1081 | case -3: |
---|
| 1082 | printf("--> Hmin/Hmax/Hstart must be positive"); |
---|
| 1083 | break; |
---|
| 1084 | case -4: |
---|
| 1085 | printf("--> Improper values for FacMin/FacMax/FacSafe/FacRej"); |
---|
| 1086 | break; |
---|
| 1087 | case -5: |
---|
| 1088 | printf("--> Improper value for ThetaMin"); |
---|
| 1089 | break; |
---|
| 1090 | case -6: |
---|
| 1091 | printf("--> Newton stopping tolerance too small"); |
---|
| 1092 | break; |
---|
| 1093 | case -7: |
---|
| 1094 | printf("--> Improper values for Qmin, Qmax"); |
---|
| 1095 | break; |
---|
| 1096 | case -8: |
---|
| 1097 | printf("--> Tolerances are too small"); |
---|
| 1098 | break; |
---|
| 1099 | case -9: |
---|
| 1100 | printf("--> No of steps exceeds maximum bound"); |
---|
| 1101 | break; |
---|
| 1102 | case -10: |
---|
| 1103 | printf("--> Step size too small: (T + 10*H = T) or H < Roundoff"); |
---|
| 1104 | break; |
---|
| 1105 | case -11: |
---|
| 1106 | printf("--> Matrix is repeatedly singular"); |
---|
| 1107 | break; |
---|
| 1108 | case -12: |
---|
| 1109 | printf("--> Non-convergence of Newton iterations"); |
---|
| 1110 | break; |
---|
| 1111 | case -13: |
---|
| 1112 | printf("--> Requested RK method not implemented"); |
---|
| 1113 | break; |
---|
| 1114 | default: |
---|
| 1115 | printf("Unknown Error code: %d \n", Code); |
---|
| 1116 | } /* end switch */ |
---|
| 1117 | |
---|
| 1118 | printf("\n T=%e, H =%e", T, H); |
---|
| 1119 | printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"); |
---|
| 1120 | |
---|
| 1121 | } /* RK_ErrorMsg */ |
---|
| 1122 | |
---|
| 1123 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1124 | /* Handles all error messages and returns SCAL */ |
---|
| 1125 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1126 | void RK_ErrorScale(int N, |
---|
| 1127 | /*~~~> Input arguments: */ |
---|
| 1128 | int ITOL, KPP_REAL AbsTol[], KPP_REAL RelTol[], KPP_REAL Y[], |
---|
| 1129 | /*~~~> Output arguments: */ |
---|
| 1130 | KPP_REAL SCAL[]) |
---|
| 1131 | { |
---|
| 1132 | /*printf("Starting RK_ErrorScale\n");*/ |
---|
| 1133 | int i; |
---|
| 1134 | if (ITOL == 0) { |
---|
| 1135 | for (i = 0; i < N; i++) { |
---|
| 1136 | SCAL[i] = ONE/(AbsTol[0]+RelTol[0]*ABS(Y[i])); |
---|
| 1137 | } /* end for loop*/ |
---|
| 1138 | } |
---|
| 1139 | else { |
---|
| 1140 | for (i = 0; i< N; i++) { |
---|
| 1141 | SCAL[i] = ONE/(AbsTol[i]+RelTol[i]*ABS(Y[i])); |
---|
| 1142 | } /* end for loop */ |
---|
| 1143 | } /* end if */ |
---|
| 1144 | /*printf("Ending RK_ErrorScale\n");*/ |
---|
| 1145 | } /* RK_ErrorScale */ |
---|
| 1146 | |
---|
| 1147 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1148 | void RK_Transform(int N, KPP_REAL Tr[][3], KPP_REAL Z1[], KPP_REAL Z2[], |
---|
| 1149 | KPP_REAL Z3[], KPP_REAL W1[], KPP_REAL W2[], KPP_REAL W3[]) |
---|
| 1150 | --> W <-- Tr x Z |
---|
| 1151 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1152 | { |
---|
| 1153 | int i; |
---|
| 1154 | KPP_REAL x1,x2,x3; |
---|
| 1155 | for (i = 0; i < N; i++) |
---|
| 1156 | { |
---|
| 1157 | x1 = Z1[i]; |
---|
| 1158 | x2 = Z2[i]; |
---|
| 1159 | x3 = Z3[i]; |
---|
| 1160 | W1[i] = Tr[0][0]*x1 + Tr[0][1]*x2 + Tr[0][2]*x3; |
---|
| 1161 | W2[i] = Tr[1][0]*x1 + Tr[1][1]*x2 + Tr[1][2]*x3; |
---|
| 1162 | W1[i] = Tr[2][0]*x1 + Tr[2][1]*x2 + Tr[2][2]*x3; |
---|
| 1163 | } end for loop |
---|
| 1164 | } end RK_Transform */ |
---|
| 1165 | |
---|
| 1166 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1167 | /*--> Constructs or evaluates a quadratic polynomial that |
---|
| 1168 | interpolates the Z solution at current step and |
---|
| 1169 | provides starting values for the next step |
---|
| 1170 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1171 | void RK_Interpolate(char action[], int N, KPP_REAL H, KPP_REAL Hold, KPP_REAL Z1[], |
---|
| 1172 | KPP_REAL Z2[], KPP_REAL Z3[], KPP_REAL CONT[][RKmax]) |
---|
| 1173 | { |
---|
| 1174 | int i; |
---|
| 1175 | KPP_REAL r,x1,x2,x3,den; |
---|
| 1176 | |
---|
| 1177 | /* Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3 */ |
---|
| 1178 | if (action == "make") { |
---|
| 1179 | den = (rkC[2]-rkC[1])*(rkC[1]-rkC[0])*(rkC[0]-rkC[2]); |
---|
| 1180 | for (i = 0; i < N; i++) |
---|
| 1181 | { |
---|
| 1182 | CONT[i][0] = (-rkC[2]*rkC[2]*rkC[1]*Z1[i] |
---|
| 1183 | +Z3[i]*rkC[1]*rkC[0]*rkC[0] |
---|
| 1184 | +rkC[1]*rkC[1]*rkC[2]*Z1[i] |
---|
| 1185 | -rkC[1]*rkC[1]*rkC[0]*Z3[i] |
---|
| 1186 | +rkC[2]*rkC[2]*rkC[0]*Z2[i] |
---|
| 1187 | -Z2[i]*rkC[2]*rkC[0]*rkC[0])/den-Z3[i]; |
---|
| 1188 | CONT[i][1] = -( rkC[0]*rkC[0]*(Z3[i]-Z2[i]) |
---|
| 1189 | + rkC[1]*rkC[1]*(Z1[i]-Z3[i]) |
---|
| 1190 | + rkC[2]*rkC[2]*(Z2[i]-Z1[i]) )/den; |
---|
| 1191 | CONT[i][2] = ( rkC[0]*(Z3[i]-Z2[i]) |
---|
| 1192 | + rkC[1]*(Z1[i]-Z3[i]) |
---|
| 1193 | + rkC[2]*(Z2[i]-Z1[i]) )/den; |
---|
| 1194 | } /* end for loop */ |
---|
| 1195 | } |
---|
| 1196 | /* Evaluate quadratic polynomial */ |
---|
| 1197 | else if (action == "eval") { |
---|
| 1198 | r = H / Hold; |
---|
| 1199 | x1 = ONE + rkC[0]*r; |
---|
| 1200 | x2 = ONE + rkC[1]*r; |
---|
| 1201 | x3 = ONE + rkC[2]*r; |
---|
| 1202 | for (i = 0; i < N; i++) |
---|
| 1203 | { |
---|
| 1204 | Z1[i] = CONT[i][0]+x1*(CONT[i][1]+x1*CONT[i][2]); |
---|
| 1205 | Z2[i] = CONT[i][0]+x2*(CONT[i][1]+x2*CONT[i][2]); |
---|
| 1206 | Z3[i] = CONT[i][0]+x3*(CONT[i][1]+x3*CONT[i][2]); |
---|
| 1207 | } /* end for loop */ |
---|
| 1208 | } /* end if */ |
---|
| 1209 | } /* RK_Interpolate */ |
---|
| 1210 | |
---|
| 1211 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1212 | /*~~~> Prepare the right-hand side for Newton iterations |
---|
| 1213 | R = Z - hA x F |
---|
| 1214 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1215 | void RK_PrepareRHS(int N, KPP_REAL T, KPP_REAL H, KPP_REAL Y[], KPP_REAL FO[], |
---|
| 1216 | KPP_REAL Z1[], KPP_REAL Z2[], KPP_REAL Z3[], KPP_REAL R1[], KPP_REAL R2[], |
---|
| 1217 | KPP_REAL R3[]) |
---|
| 1218 | { |
---|
| 1219 | KPP_REAL TMP[N], F[N]; |
---|
| 1220 | WCOPY(N,Z1,1,R1,1); /* R1 <- Z1 */ |
---|
| 1221 | WCOPY(N,Z2,1,R2,1); /* R2 <- Z2 */ |
---|
| 1222 | WCOPY(N,Z3,1,R3,1); /* R3 <- Z3 */ |
---|
| 1223 | |
---|
| 1224 | if (rkMethod == L3A) |
---|
| 1225 | { |
---|
| 1226 | WAXPY(N,-H*rkA[0][0],FO,1,R1,1); /* R1 <- R1 - h*A_10*FO */ |
---|
| 1227 | WAXPY(N,-H*rkA[1][0],FO,1,R2,1); /* R2 <- R2 - h*A_20*FO */ |
---|
| 1228 | WAXPY(N,-H*rkA[2][0],FO,1,R3,1); /* R3 <- R3 - h*A_30*FO */ |
---|
| 1229 | } /* end if */ |
---|
| 1230 | |
---|
| 1231 | WADD(N,Y,Z1,TMP); /* TMP <- Y + Z1 */ |
---|
| 1232 | FUN_CHEM(T+rkC[0]*H,TMP,F); /* F1 <- Fun(Y+Z1) */ |
---|
| 1233 | WAXPY(N,-H*rkA[0][0],F,1,R1,1); /* R1 <- R1 - h*A_11*F1 */ |
---|
| 1234 | WAXPY(N,-H*rkA[1][0],F,1,R2,1); /* R2 <- R2 - h*A_21*F1 */ |
---|
| 1235 | WAXPY(N,-H*rkA[2][0],F,1,R3,1); /* R3 <- R3 - h*A_31*F1 */ |
---|
| 1236 | |
---|
| 1237 | WADD(N,Y,Z2,TMP); /* TMP <- Y + Z2 */ |
---|
| 1238 | FUN_CHEM(T+rkC[1]*H,TMP,F); /* F2 <- Fun(Y+Z2) */ |
---|
| 1239 | WAXPY(N,-H*rkA[0][1],F,1,R1,1); /* R1 <- R1 - h*A_12*F2 */ |
---|
| 1240 | WAXPY(N,-H*rkA[1][1],F,1,R2,1); /* R2 <- R2 - h*A_22*F2 */ |
---|
| 1241 | WAXPY(N,-H*rkA[2][1],F,1,R3,1); /* R3 <- R3 - h*A_32*F2 */ |
---|
| 1242 | |
---|
| 1243 | WADD(N,Y,Z3,TMP); /* TMP <- Y + Z3 */ |
---|
| 1244 | FUN_CHEM(T+rkC[2]*H,TMP,F); /* F3 <- Fun(Y+Z3) */ |
---|
| 1245 | WAXPY(N,-H*rkA[0][2],F,1,R1,1); /* R1 <- R1 - h*A_13*F3 */ |
---|
| 1246 | WAXPY(N,-H*rkA[1][2],F,1,R2,1); /* R2 <- R2 - h*A_23*F3 */ |
---|
| 1247 | WAXPY(N,-H*rkA[2][2],F,1,R3,1); /* R3 <- R3 - h*A_33*F3 */ |
---|
| 1248 | } /* RK_PrepareRHS */ |
---|
| 1249 | |
---|
| 1250 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1251 | /*~~~> Compute the matrices E1 and E2 and their decompositions |
---|
| 1252 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1253 | void RK_Decomp(int N, KPP_REAL H, KPP_REAL FJAC[], KPP_REAL E1[], |
---|
| 1254 | int IP1[], KPP_REAL E2R[], KPP_REAL E2I[], int IP2[], |
---|
| 1255 | int* ISING, int ISTATUS[]) |
---|
| 1256 | { |
---|
| 1257 | /*printf("Starting RK_Decomp\n");*/ |
---|
| 1258 | KPP_REAL Alpha, Beta, Gamma; |
---|
| 1259 | int i,j; |
---|
| 1260 | |
---|
| 1261 | Gamma = rkGamma / H; |
---|
| 1262 | Alpha = rkAlpha / H; |
---|
| 1263 | Beta = rkBeta / H; |
---|
| 1264 | |
---|
| 1265 | for (i = 0; i < LU_NONZERO; i++) { |
---|
| 1266 | E1[i] = -FJAC[i]; |
---|
| 1267 | } |
---|
| 1268 | for (i = 0; i < NVAR; i++) { |
---|
| 1269 | j = LU_DIAG[i]; |
---|
| 1270 | E1[j] = E1[j] + Gamma; |
---|
| 1271 | } |
---|
| 1272 | *ISING = KppDecomp(E1); |
---|
| 1273 | /*~~~> Note: for a full matrix use Lapack: |
---|
| 1274 | for (j = 0; j < N; j++) { |
---|
| 1275 | for (i = 0; i < N; i++) |
---|
| 1276 | { |
---|
| 1277 | E1[i,j] = -FJAC[i][j]; |
---|
| 1278 | } |
---|
| 1279 | E1[i][j] = E1[i][j]+Gamma; |
---|
| 1280 | } |
---|
| 1281 | DGETRF(N, N, E1, N, IP1, ISING); */ |
---|
| 1282 | |
---|
| 1283 | if (*ISING != 0) { |
---|
| 1284 | /*printf("Matrix is singular, ISING=%d\n", *ISING);*/ |
---|
| 1285 | ISTATUS[Ndec]++; |
---|
| 1286 | return; |
---|
| 1287 | } |
---|
| 1288 | |
---|
| 1289 | for (i = 0; i < LU_NONZERO; i++) { |
---|
| 1290 | E2R[i] = (KPP_REAL)(-FJAC[i]); |
---|
| 1291 | E2I[i] = ZERO; |
---|
| 1292 | } |
---|
| 1293 | for (i = 0; i < NVAR; i++) { |
---|
| 1294 | j = LU_DIAG[i]; |
---|
| 1295 | E2R[j] = E2R[j] + Alpha; |
---|
| 1296 | E2I[j] = E2I[j] + Beta; |
---|
| 1297 | } |
---|
| 1298 | *ISING = KppDecompCmplxR(E2R, E2I); |
---|
| 1299 | /*printf("Matrix is singular, ISING=%d\n", *ISING);*/ |
---|
| 1300 | /*~~~> Note: for a full matrix use Lapack: |
---|
| 1301 | for (j = 0; j < N; j++) { |
---|
| 1302 | for (i = 0; i < N; i++) { |
---|
| 1303 | E2[i,j] = DCMPLX( -FJAC[i,j], ZERO ); |
---|
| 1304 | } |
---|
| 1305 | E2[j,j] = E2[j,j] + CMPLX( Alpha, Beta ); |
---|
| 1306 | } |
---|
| 1307 | ZGETRF(N,N,E2,N,IP2,ISING); */ |
---|
| 1308 | ISTATUS[Ndec]++; |
---|
| 1309 | } /*RK_Decomp */ |
---|
| 1310 | |
---|
| 1311 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1312 | void RK_Solve(int N, KPP_REAL H, KPP_REAL E1[], int IP1[], KPP_REAL E2R[], |
---|
| 1313 | KPP_REAL E2I[], int IP2[],KPP_REAL R1[], KPP_REAL R2[], KPP_REAL R3[], |
---|
| 1314 | int ISTATUS[]) |
---|
| 1315 | { |
---|
| 1316 | KPP_REAL x1, x2, x3; |
---|
| 1317 | KPP_REAL BCR[N], BCI[N]; |
---|
| 1318 | int i; |
---|
| 1319 | |
---|
| 1320 | /* Z <- h^{-1) T^{-1) A^{-1) x Z */ |
---|
| 1321 | for (i = 0; i < N; i++) |
---|
| 1322 | { |
---|
| 1323 | x1 = R1[i]/H; |
---|
| 1324 | x2 = R2[i]/H; |
---|
| 1325 | x3 = R3[i]/H; |
---|
| 1326 | R1[i] = rkTinvAinv[0][0]*x1 + rkTinvAinv[0][1]*x2 + rkTinvAinv[0][2]*x3; |
---|
| 1327 | R2[i] = rkTinvAinv[1][0]*x1 + rkTinvAinv[1][1]*x2 + rkTinvAinv[1][2]*x3; |
---|
| 1328 | R3[i] = rkTinvAinv[2][0]*x1 + rkTinvAinv[2][1]*x2 + rkTinvAinv[2][2]*x3; |
---|
| 1329 | } |
---|
| 1330 | KppSolve(E1,R1); |
---|
| 1331 | /*~~~> Note: for a full matrix use Lapack: |
---|
| 1332 | DGETRS('N',5,1,E1,N,IP1,R1,5,0); */ |
---|
| 1333 | for(i = 0; i < N; i++) |
---|
| 1334 | { |
---|
| 1335 | BCR[i] = (KPP_REAL)(R2[i]); |
---|
| 1336 | BCI[i] = (KPP_REAL)(R3[i]); |
---|
| 1337 | } |
---|
| 1338 | KppSolveCmplxR(E2R,E2I,BCR,BCI); |
---|
| 1339 | /*~~~> Note: for a full matrix use Lapack: |
---|
| 1340 | ZGETRS ('N',N,1,E2,N,IP2,BC,N,0); */ |
---|
| 1341 | for (i = 0; i < N; i++) |
---|
| 1342 | { |
---|
| 1343 | R2[i] = (KPP_REAL)( BCR[i] ); |
---|
| 1344 | R3[i] = (KPP_REAL)( BCI[i] ); |
---|
| 1345 | } |
---|
| 1346 | |
---|
| 1347 | /* Z <- T x Z */ |
---|
| 1348 | for (i = 0; i < N; i++) |
---|
| 1349 | { |
---|
| 1350 | x1 = R1[i]; |
---|
| 1351 | x2 = R2[i]; |
---|
| 1352 | x3 = R3[i]; |
---|
| 1353 | R1[i] = rkT[0][0]*x1 + rkT[0][1]*x2 + rkT[0][2]*x3; |
---|
| 1354 | R2[i] = rkT[1][0]*x1 + rkT[1][1]*x2 + rkT[1][2]*x3; |
---|
| 1355 | R3[i] = rkT[2][0]*x1 + rkT[2][1]*x2 + rkT[2][2]*x3; |
---|
| 1356 | } |
---|
| 1357 | ISTATUS[Nsol]++; |
---|
| 1358 | } /* RK_Solve */ |
---|
| 1359 | |
---|
| 1360 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1361 | void RK_ErrorEstimate(int N, KPP_REAL H, KPP_REAL T, KPP_REAL Y[], |
---|
| 1362 | KPP_REAL FO[], KPP_REAL E1[], int IP1[], KPP_REAL Z1[], |
---|
| 1363 | KPP_REAL Z2[], KPP_REAL Z3[], KPP_REAL SCAL[], KPP_REAL* Err, |
---|
| 1364 | int FirstStep, int Reject, int ISTATUS[]) |
---|
| 1365 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1366 | { |
---|
| 1367 | KPP_REAL F1[N],F2[N],TMP[N]; |
---|
| 1368 | KPP_REAL HrkE1,HrkE2,HrkE3; |
---|
| 1369 | int i; |
---|
| 1370 | |
---|
| 1371 | HrkE1 = rkE[0]/H; |
---|
| 1372 | HrkE2 = rkE[1]/H; |
---|
| 1373 | HrkE3 = rkE[2]/H; |
---|
| 1374 | |
---|
| 1375 | for (i = 0; i < N; i++) |
---|
| 1376 | { |
---|
| 1377 | F2[i] = HrkE1*Z1[i]+HrkE2*Z2[i]+HrkE3*Z3[i]; |
---|
| 1378 | TMP[i] = rkE[0]*FO[i] + F2[i]; |
---|
| 1379 | } |
---|
| 1380 | |
---|
| 1381 | KppSolve(E1, TMP); |
---|
| 1382 | if ((rkMethod == R1A) || (rkMethod == GAU) || (rkMethod == L3A)) |
---|
| 1383 | KppSolve(E1,TMP); |
---|
| 1384 | if (rkMethod == GAU) |
---|
| 1385 | KppSolve(E1,TMP); |
---|
| 1386 | /*~~~> Note: for a full matrix use Lapack: |
---|
| 1387 | DGETRS ('N',N,1,E1,N,IP1,TMP,N,0); |
---|
| 1388 | if ((rkMethod == R1A) || (rkMethod == GAU) || (rkMethod == L3A)) |
---|
| 1389 | DGETRS('N',N,1,E1,N,IP1,TMP,N,0); |
---|
| 1390 | if (rkMethod == GAU) |
---|
| 1391 | DGETRS('N',N,1,E1,N,IP1,TMP,N,0); */ |
---|
| 1392 | |
---|
| 1393 | *Err = RK_ErrorNorm(N,SCAL,TMP); |
---|
| 1394 | |
---|
| 1395 | if (*Err < ONE) |
---|
| 1396 | return; |
---|
| 1397 | /* firej */ |
---|
| 1398 | if ((FirstStep == 1) || (Reject == 1)) { |
---|
| 1399 | for (i = 0; i < N; i++) { |
---|
| 1400 | TMP[i] = Y[i] + TMP[i]; |
---|
| 1401 | } /* end for */ |
---|
| 1402 | FUN_CHEM(T,TMP,F1); |
---|
| 1403 | ISTATUS[Nfun]++; |
---|
| 1404 | for (i = 0; i < N; i++) { |
---|
| 1405 | TMP[i] = F1[i] + F2[i]; |
---|
| 1406 | } /* end for */ |
---|
| 1407 | |
---|
| 1408 | KppSolve(E1, TMP); |
---|
| 1409 | /*~~~> Note: for a full matrix use Lapack: |
---|
| 1410 | DGETRS ('N',N,1,E1,N,IP1,TMP,N,0); */ |
---|
| 1411 | *Err = RK_ErrorNorm(N,SCAL,TMP); |
---|
| 1412 | } /* end if firej */ |
---|
| 1413 | } /* RK_ErrorEstimate */ |
---|
| 1414 | |
---|
| 1415 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1416 | KPP_REAL RK_ErrorNorm(int N, KPP_REAL SCAL[], KPP_REAL DY[]) |
---|
| 1417 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1418 | { |
---|
| 1419 | int i; |
---|
| 1420 | KPP_REAL RK_ErrorNorm = ZERO; |
---|
| 1421 | for (i = 0; i < N; i++) { |
---|
| 1422 | RK_ErrorNorm = RK_ErrorNorm + (DY[i]*SCAL[i]) * (DY[i]*SCAL[i]); |
---|
| 1423 | } |
---|
| 1424 | RK_ErrorNorm = MAX( SQRT(RK_ErrorNorm/N), (KPP_REAL)1.0e-10 ); |
---|
| 1425 | return RK_ErrorNorm; |
---|
| 1426 | } /* RK_ErrorNorm */ |
---|
| 1427 | |
---|
| 1428 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1429 | The coefficients of the 3-stage Radau-2A method |
---|
| 1430 | (given to ~30 accurate digits) |
---|
| 1431 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1432 | void Radau2A_Coefficients() |
---|
| 1433 | { |
---|
| 1434 | /*printf("Starting Radau2A\n");*/ |
---|
| 1435 | /* The coefficients of the Radau2A method */ |
---|
| 1436 | KPP_REAL b0; |
---|
| 1437 | /* b0 = (KPP_REAL)1.0; */ |
---|
| 1438 | if (SdirkError == 1) { |
---|
| 1439 | b0 = (KPP_REAL)(0.2e-01); |
---|
| 1440 | /*printf("b0 = %f\n", b0);*/ |
---|
| 1441 | } |
---|
| 1442 | else { |
---|
| 1443 | b0 = (KPP_REAL)(0.5e-01); |
---|
| 1444 | /*printf("b0 = %f\n", b0);*/ |
---|
| 1445 | } /*end if */ |
---|
| 1446 | /* The coefficients of the Radau2A method */ |
---|
| 1447 | rkMethod = R2A; |
---|
| 1448 | |
---|
| 1449 | rkA[0][0] = (KPP_REAL)(1.968154772236604258683861429918299e-01); |
---|
| 1450 | rkA[0][1] = (KPP_REAL)(-6.55354258501983881085227825696087e-02); |
---|
| 1451 | rkA[0][2] = (KPP_REAL)(2.377097434822015242040823210718965e-02); |
---|
| 1452 | rkA[1][0] = (KPP_REAL)(3.944243147390872769974116714584975e-01); |
---|
| 1453 | rkA[1][1] = (KPP_REAL)(2.920734116652284630205027458970589e-01); |
---|
| 1454 | rkA[1][2] = (KPP_REAL)(-4.154875212599793019818600988496743e-02); |
---|
| 1455 | rkA[2][0] = (KPP_REAL)(3.764030627004672750500754423692808e-01); |
---|
| 1456 | rkA[2][1] = (KPP_REAL)(5.124858261884216138388134465196080e-01); |
---|
| 1457 | rkA[2][2] = (KPP_REAL)(1.111111111111111111111111111111111e-01); |
---|
| 1458 | |
---|
| 1459 | rkB[0] = (KPP_REAL)(3.764030627004672750500754423692808e-01); |
---|
| 1460 | rkB[1] = (KPP_REAL)(5.124858261884216138388134465196080e-01); |
---|
| 1461 | rkB[2] = (KPP_REAL)(1111111111111111111111111111111111e-01); |
---|
| 1462 | |
---|
| 1463 | rkC[0] = (KPP_REAL)(1.550510257216821901802715925294109e-01); |
---|
| 1464 | rkC[1] = (KPP_REAL)(6.449489742783178098197284074705891e-01); |
---|
| 1465 | rkC[2] = ONE; |
---|
| 1466 | |
---|
| 1467 | /* New solution: H* Sum B_j*f(Z_j) = Sum D_j*Z_j */ |
---|
| 1468 | rkD[0] = ZERO; |
---|
| 1469 | rkD[1] = ZERO; |
---|
| 1470 | rkD[2] = ONE; |
---|
| 1471 | |
---|
| 1472 | /* Classical error estimator: */ |
---|
| 1473 | /* H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j */ |
---|
| 1474 | rkE[0] = ONE*b0; |
---|
| 1475 | rkE[1] = (KPP_REAL)(-10.04880939982741556246032950764708)*b0; |
---|
| 1476 | rkE[2] = (KPP_REAL)1.382142733160748895793662840980412*b0; |
---|
| 1477 | rkE[3] = (KPP_REAL)(-.3333333333333333333333333333333333)*b0; |
---|
| 1478 | |
---|
| 1479 | /* Sdirk error estimator */ |
---|
| 1480 | rkBgam[0] = b0; |
---|
| 1481 | rkBgam[1] = (KPP_REAL).3764030627004672750500754423692807 |
---|
| 1482 | -(KPP_REAL)1.558078204724922382431975370686279*b0; |
---|
| 1483 | rkBgam[2] = (KPP_REAL).8914115380582557157653087040196118*b0 |
---|
| 1484 | +(KPP_REAL).5124858261884216138388134465196077; |
---|
| 1485 | rkBgam[3] = (KPP_REAL)(-.1637777184845662566367174924883037) |
---|
| 1486 | -(KPP_REAL).3333333333333333333333333333333333*b0; |
---|
| 1487 | rkBgam[4] = (KPP_REAL).2748888295956773677478286035994148; |
---|
| 1488 | |
---|
| 1489 | /* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */ |
---|
| 1490 | rkTheta[0] = (KPP_REAL)(-1.520677486405081647234271944611547) |
---|
| 1491 | -(KPP_REAL)10.04880939982741556246032950764708*b0; |
---|
| 1492 | rkTheta[1] = (KPP_REAL)2.070455145596436382729929151810376 |
---|
| 1493 | +(KPP_REAL)1.382142733160748895793662840980413*b0; |
---|
| 1494 | rkTheta[2] = (KPP_REAL)(-.3333333333333333333333333333333333)*b0 |
---|
| 1495 | -(KPP_REAL).3744441479783868387391430179970741; |
---|
| 1496 | |
---|
| 1497 | /* Local order of error estimator */ |
---|
| 1498 | if ( b0 == ZERO) |
---|
| 1499 | rkELO = (KPP_REAL)6.0; |
---|
| 1500 | else |
---|
| 1501 | rkELO = (KPP_REAL)4.0; |
---|
| 1502 | |
---|
| 1503 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1504 | ~~~> Diagonalize the RK matrix: |
---|
| 1505 | rkTinv * inv(rkA) * rkT = |
---|
| 1506 | | rkGamma 0 0 | |
---|
| 1507 | | 0 rkAlpha -rkBeta | |
---|
| 1508 | | 0 rkBeta rkAlpha | |
---|
| 1509 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1510 | |
---|
| 1511 | rkGamma = (KPP_REAL)3.637834252744495732208418513577775; |
---|
| 1512 | rkAlpha = (KPP_REAL)2.681082873627752133895790743211112; |
---|
| 1513 | rkBeta = (KPP_REAL)3.050430199247410569426377624787569; |
---|
| 1514 | |
---|
| 1515 | rkT[0][0] = (KPP_REAL)(9.443876248897524148749007950641664e-02); |
---|
| 1516 | rkT[0][1] = (KPP_REAL)(-1.412552950209542084279903838077973e-01); |
---|
| 1517 | rkT[0][2] = (KPP_REAL)(-3.00291941051474244918611170890539e-02); |
---|
| 1518 | rkT[1][0] = (KPP_REAL)(2.502131229653333113765090675125018e-01); |
---|
| 1519 | rkT[1][1] = (KPP_REAL)(2.041293522937999319959908102983381e-01); |
---|
| 1520 | rkT[1][2] = (KPP_REAL)(3.829421127572619377954382335998733e-01); |
---|
| 1521 | rkT[2][0] = ONE; |
---|
| 1522 | rkT[2][1] = ONE; |
---|
| 1523 | rkT[2][2] = ZERO; |
---|
| 1524 | |
---|
| 1525 | rkTinv[0][0] = (KPP_REAL)4.178718591551904727346462658512057; |
---|
| 1526 | rkTinv[0][1] = (KPP_REAL)(3.27682820761062387082533272429617e-01); |
---|
| 1527 | rkTinv[0][2] = (KPP_REAL)(5.233764454994495480399309159089876e-01); |
---|
| 1528 | rkTinv[1][0] = (KPP_REAL)(-4.178718591551904727346462658512057); |
---|
| 1529 | rkTinv[1][1] = (KPP_REAL)(-3.27682820761062387082533272429617e-01); |
---|
| 1530 | rkTinv[1][2] = (KPP_REAL)(4.766235545005504519600690840910124e-01); |
---|
| 1531 | rkTinv[2][0] = (KPP_REAL)(-5.02872634945786875951247343139544e-01); |
---|
| 1532 | rkTinv[2][1] = (KPP_REAL)2.571926949855605429186785353601676; |
---|
| 1533 | rkTinv[2][2] = (KPP_REAL)(-5.960392048282249249688219110993024e-01); |
---|
| 1534 | |
---|
| 1535 | rkTinvAinv[0][0] = (KPP_REAL)(1.520148562492775501049204957366528e+01); |
---|
| 1536 | rkTinvAinv[0][1] = (KPP_REAL)1.192055789400527921212348994770778; |
---|
| 1537 | rkTinvAinv[0][2] = (KPP_REAL)1.903956760517560343018332287285119; |
---|
| 1538 | rkTinvAinv[1][0] = (KPP_REAL)(-9.669512977505946748632625374449567); |
---|
| 1539 | rkTinvAinv[1][1] = (KPP_REAL)(-8.724028436822336183071773193986487); |
---|
| 1540 | rkTinvAinv[1][2] = (KPP_REAL)3.096043239482439656981667712714881; |
---|
| 1541 | rkTinvAinv[2][0] = (KPP_REAL)(-1.409513259499574544876303981551774e+01); |
---|
| 1542 | rkTinvAinv[2][1] = (KPP_REAL)5.895975725255405108079130152868952; |
---|
| 1543 | rkTinvAinv[2][2] = (KPP_REAL)(-1.441236197545344702389881889085515e-01); |
---|
| 1544 | |
---|
| 1545 | rkAinvT[0][0] = (KPP_REAL).3435525649691961614912493915818282; |
---|
| 1546 | rkAinvT[0][1] = (KPP_REAL)(-.4703191128473198422370558694426832); |
---|
| 1547 | rkAinvT[0][2] = (KPP_REAL).3503786597113668965366406634269080; |
---|
| 1548 | rkAinvT[1][0] = (KPP_REAL).9102338692094599309122768354288852; |
---|
| 1549 | rkAinvT[1][1] = (KPP_REAL)1.715425895757991796035292755937326; |
---|
| 1550 | rkAinvT[1][2] = (KPP_REAL).4040171993145015239277111187301784; |
---|
| 1551 | rkAinvT[2][0] = (KPP_REAL)3.637834252744495732208418513577775; |
---|
| 1552 | rkAinvT[2][1] = (KPP_REAL)2.681082873627752133895790743211112; |
---|
| 1553 | rkAinvT[2][2] = (KPP_REAL)(-3.050430199247410569426377624787569); |
---|
| 1554 | /*printf("Ending Radau2A");*/ |
---|
| 1555 | } /* Radau2A_Coefficients */ |
---|
| 1556 | |
---|
| 1557 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1558 | The coefficients of the 3-stage Lobatto-3C method |
---|
| 1559 | (given to ~30 accurate digits) |
---|
| 1560 | The parameter b0 can be chosen to tune the error estimator |
---|
| 1561 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1562 | void Lobatto3C_Coefficients () |
---|
| 1563 | { |
---|
| 1564 | KPP_REAL b0; |
---|
| 1565 | rkMethod = L3C; |
---|
| 1566 | /* b0 = 1.0d0 */ |
---|
| 1567 | if (SdirkError == 1) |
---|
| 1568 | b0 = (KPP_REAL)0.2; |
---|
| 1569 | else |
---|
| 1570 | b0 = (KPP_REAL)0.5; |
---|
| 1571 | |
---|
| 1572 | /* The coefficients of the Lobatto3C method */ |
---|
| 1573 | rkA[0][0] = (KPP_REAL).1666666666666666666666666666666667; |
---|
| 1574 | rkA[0][1] = (KPP_REAL)(-.3333333333333333333333333333333333); |
---|
| 1575 | rkA[0][2] = (KPP_REAL).1666666666666666666666666666666667; |
---|
| 1576 | rkA[1][0] = (KPP_REAL).1666666666666666666666666666666667; |
---|
| 1577 | rkA[1][1] = (KPP_REAL).4166666666666666666666666666666667; |
---|
| 1578 | rkA[1][2] = (KPP_REAL)(-.8333333333333333333333333333333333e-01); |
---|
| 1579 | rkA[2][0] = (KPP_REAL).1666666666666666666666666666666667; |
---|
| 1580 | rkA[2][1] = (KPP_REAL).6666666666666666666666666666666667; |
---|
| 1581 | rkA[2][2] = (KPP_REAL).1666666666666666666666666666666667; |
---|
| 1582 | |
---|
| 1583 | rkB[0] = (KPP_REAL).1666666666666666666666666666666667; |
---|
| 1584 | rkB[1] = (KPP_REAL).6666666666666666666666666666666667; |
---|
| 1585 | rkB[2] = (KPP_REAL).1666666666666666666666666666666667; |
---|
| 1586 | |
---|
| 1587 | rkC[0] = ZERO; |
---|
| 1588 | rkC[1] = (KPP_REAL)0.5; |
---|
| 1589 | rkC[2] = ONE; |
---|
| 1590 | |
---|
| 1591 | /* Classical error estimator, embedded solution: */ |
---|
| 1592 | rkBhat[0] = b0; |
---|
| 1593 | rkBhat[1] = (KPP_REAL).16666666666666666666666666666666667-b0; |
---|
| 1594 | rkBhat[2] = (KPP_REAL).66666666666666666666666666666666667; |
---|
| 1595 | rkBhat[3] = (KPP_REAL).16666666666666666666666666666666667; |
---|
| 1596 | |
---|
| 1597 | /* New solution: h Sum_j b_j f(Z_j) = sum d_j Z_j */ |
---|
| 1598 | rkD[0] = ZERO; |
---|
| 1599 | rkD[1] = ZERO; |
---|
| 1600 | rkD[2] = ONE; |
---|
| 1601 | |
---|
| 1602 | /* Classical error estimator: */ |
---|
| 1603 | /* H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j */ |
---|
| 1604 | rkE[0] = (KPP_REAL).3808338772072650364017425226487022*b0; |
---|
| 1605 | rkE[1] = (KPP_REAL)(-1.142501631621795109205227567946107)*b0; |
---|
| 1606 | rkE[2] = (KPP_REAL)(-1.523335508829060145606970090594809)*b0; |
---|
| 1607 | rkE[3] = (KPP_REAL).3808338772072650364017425226487022*b0; |
---|
| 1608 | |
---|
| 1609 | /* Sdirk error estimator */ |
---|
| 1610 | rkBgam[0] = b0; |
---|
| 1611 | rkBgam[1] = (KPP_REAL).1666666666666666666666666666666667-1.0*b0; |
---|
| 1612 | rkBgam[2] = (KPP_REAL).6666666666666666666666666666666667; |
---|
| 1613 | rkBgam[3] = (KPP_REAL)(-.2141672105405983697350758559820354); |
---|
| 1614 | rkBgam[4] = (KPP_REAL).3808338772072650364017425226487021; |
---|
| 1615 | |
---|
| 1616 | /* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */ |
---|
| 1617 | rkTheta[0] = -3.0*b0-(KPP_REAL).3808338772072650364017425226487021; |
---|
| 1618 | rkTheta[1] = -4.0*b0+(KPP_REAL)1.523335508829060145606970090594808; |
---|
| 1619 | rkTheta[2] = (KPP_REAL)(-.142501631621795109205227567946106)+b0; |
---|
| 1620 | |
---|
| 1621 | /* Local order of error estimator */ |
---|
| 1622 | if (b0 == ZERO) |
---|
| 1623 | rkELO = 5.0; |
---|
| 1624 | else |
---|
| 1625 | rkELO = 4.0; |
---|
| 1626 | |
---|
| 1627 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1628 | ~~~> Diagonalize the RK matrix: |
---|
| 1629 | rkTinv * inv(rkA) * rkT = |
---|
| 1630 | | rkGamma 0 0 | |
---|
| 1631 | | 0 rkAlpha -rkBeta | |
---|
| 1632 | | 0 rkBeta rkAlpha | |
---|
| 1633 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1634 | |
---|
| 1635 | rkGamma = (KPP_REAL)2.625816818958466716011888933765284; |
---|
| 1636 | rkAlpha = (KPP_REAL)1.687091590520766641994055533117359; |
---|
| 1637 | rkBeta = (KPP_REAL)2.508731754924880510838743672432351; |
---|
| 1638 | |
---|
| 1639 | rkT[0][1] = ONE; |
---|
| 1640 | rkT[0][1] = ONE; |
---|
| 1641 | rkT[0][2] = ZERO; |
---|
| 1642 | rkT[1][0] = (KPP_REAL).4554100411010284672111720348287483; |
---|
| 1643 | rkT[1][1] = (KPP_REAL)(-.6027050205505142336055860174143743); |
---|
| 1644 | rkT[1][2] = (KPP_REAL)(-.4309321229203225731070721341350346); |
---|
| 1645 | rkT[2][0] = (KPP_REAL)2.195823345445647152832799205549709; |
---|
| 1646 | rkT[2][1] = (KPP_REAL)(-1.097911672722823576416399602774855); |
---|
| 1647 | rkT[2][2] = (KPP_REAL).7850032632435902184104551358922130; |
---|
| 1648 | |
---|
| 1649 | rkTinv[0][0] = (KPP_REAL).4205559181381766909344950150991349; |
---|
| 1650 | rkTinv[0][1] = (KPP_REAL).3488903392193734304046467270632057; |
---|
| 1651 | rkTinv[0][2] = (KPP_REAL).1915253879645878102698098373933487; |
---|
| 1652 | rkTinv[1][0] = (KPP_REAL).5794440818618233090655049849008650; |
---|
| 1653 | rkTinv[1][1] = (KPP_REAL)(-.3488903392193734304046467270632057); |
---|
| 1654 | rkTinv[1][2] = (KPP_REAL)(-.1915253879645878102698098373933487); |
---|
| 1655 | rkTinv[2][0] = (KPP_REAL)(-.3659705575742745254721332009249516); |
---|
| 1656 | rkTinv[2][1] = (KPP_REAL)(-1.463882230297098101888532803699806); |
---|
| 1657 | rkTinv[2][2] = (KPP_REAL).4702733607340189781407813565524989; |
---|
| 1658 | |
---|
| 1659 | rkTinvAinv[0][0] = (KPP_REAL)1.104302803159744452668648155627548; |
---|
| 1660 | rkTinvAinv[0][1] = (KPP_REAL).916122120694355522658740710823143; |
---|
| 1661 | rkTinvAinv[0][2] = (KPP_REAL).5029105849749601702795812241441172; |
---|
| 1662 | rkTinvAinv[1][0] = (KPP_REAL)1.895697196840255547331351844372453; |
---|
| 1663 | rkTinvAinv[1][1] = (KPP_REAL)3.083877879305644477341259289176857; |
---|
| 1664 | rkTinvAinv[1][2] = (KPP_REAL)(-1.502910584974960170279581224144117); |
---|
| 1665 | rkTinvAinv[2][0] = (KPP_REAL).8362439183082935036129145574774502; |
---|
| 1666 | rkTinvAinv[2][1] = (KPP_REAL)(-3.344975673233174014451658229909802); |
---|
| 1667 | rkTinvAinv[2][2] = (KPP_REAL).312908409479233358005944466882642; |
---|
| 1668 | |
---|
| 1669 | rkAinvT[0][0] = (KPP_REAL)2.625816818958466716011888933765282; |
---|
| 1670 | rkAinvT[0][1] = (KPP_REAL)1.687091590520766641994055533117358; |
---|
| 1671 | rkAinvT[0][2] = (KPP_REAL)(-2.508731754924880510838743672432351); |
---|
| 1672 | rkAinvT[1][0] = (KPP_REAL)1.195823345445647152832799205549710; |
---|
| 1673 | rkAinvT[1][1] = (KPP_REAL)(-2.097911672722823576416399602774855); |
---|
| 1674 | rkAinvT[1][2] = (KPP_REAL).7850032632435902184104551358922130; |
---|
| 1675 | rkAinvT[2][0] = (KPP_REAL)5.765829871932827589653709477334136; |
---|
| 1676 | rkAinvT[2][1] = (KPP_REAL).1170850640335862051731452613329320; |
---|
| 1677 | rkAinvT[2][2] = (KPP_REAL)4.078738281412060947659653944216779; |
---|
| 1678 | } /* Lobatto3C_Coefficients */ |
---|
| 1679 | |
---|
| 1680 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1681 | The coefficients of the 3-stage Gauss method |
---|
| 1682 | (given to ~30 accurate digits) |
---|
| 1683 | The parameter b3 can be chosen by the user |
---|
| 1684 | to tune the error estimator |
---|
| 1685 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1686 | void Gauss_Coefficients() |
---|
| 1687 | { |
---|
| 1688 | rkMethod = GAU; |
---|
| 1689 | |
---|
| 1690 | /* The coefficients of the Gauss method */ |
---|
| 1691 | KPP_REAL b0; |
---|
| 1692 | /*b0 = (KPP_REAL)4.0; */ |
---|
| 1693 | b0 = (KPP_REAL)0.1; |
---|
| 1694 | |
---|
| 1695 | /* The coefficients of the Gauss method */ |
---|
| 1696 | rkA[0][0] = (KPP_REAL).1388888888888888888888888888888889; |
---|
| 1697 | rkA[0][1] = (KPP_REAL)(-.359766675249389034563954710966045e-01); |
---|
| 1698 | rkA[0][2] = (KPP_REAL)(.97894440153083260495800422294756e-02); |
---|
| 1699 | rkA[1][0] = (KPP_REAL).3002631949808645924380249472131556; |
---|
| 1700 | rkA[1][1] = (KPP_REAL).2222222222222222222222222222222222; |
---|
| 1701 | rkA[1][2] = (KPP_REAL)(-.224854172030868146602471694353778e-01); |
---|
| 1702 | rkA[2][0] = (KPP_REAL).2679883337624694517281977355483022; |
---|
| 1703 | rkA[2][1] = (KPP_REAL).4804211119693833479008399155410489; |
---|
| 1704 | rkA[2][2] = (KPP_REAL).1388888888888888888888888888888889; |
---|
| 1705 | |
---|
| 1706 | rkB[0] = (KPP_REAL).2777777777777777777777777777777778; |
---|
| 1707 | rkB[1] = (KPP_REAL).4444444444444444444444444444444444; |
---|
| 1708 | rkB[2] = (KPP_REAL).2777777777777777777777777777777778; |
---|
| 1709 | |
---|
| 1710 | rkC[0] = (KPP_REAL).1127016653792583114820734600217600; |
---|
| 1711 | rkC[1] = (KPP_REAL).5000000000000000000000000000000000; |
---|
| 1712 | rkC[2] = (KPP_REAL).8872983346207416885179265399782400; |
---|
| 1713 | |
---|
| 1714 | /* Classical error estimator, embedded solution: */ |
---|
| 1715 | rkBhat[0] = b0; |
---|
| 1716 | rkBhat[1] = (KPP_REAL)(-1.4788305577012361475298775666303999)*b0 |
---|
| 1717 | +(KPP_REAL).27777777777777777777777777777777778; |
---|
| 1718 | rkBhat[2] = (KPP_REAL).44444444444444444444444444444444444 |
---|
| 1719 | +(KPP_REAL).66666666666666666666666666666666667*b0; |
---|
| 1720 | rkBhat[3] = (KPP_REAL)(-.18783610896543051913678910003626672)*b0 |
---|
| 1721 | +(KPP_REAL).27777777777777777777777777777777778; |
---|
| 1722 | |
---|
| 1723 | /* New solution: h Sum_j b_j f(Z_j) = sum d_j Z_j */ |
---|
| 1724 | rkD[0] = (KPP_REAL)(.1666666666666666666666666666666667e+01); |
---|
| 1725 | rkD[1] = (KPP_REAL)(-.1333333333333333333333333333333333e+01); |
---|
| 1726 | rkD[2] = (KPP_REAL)(.1666666666666666666666666666666667e+01); |
---|
| 1727 | |
---|
| 1728 | /* Classical error estimator: */ |
---|
| 1729 | /* H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j */ |
---|
| 1730 | rkE[0] = (KPP_REAL).2153144231161121782447335303806954*b0; |
---|
| 1731 | rkE[1] = (KPP_REAL)(-2.825278112319014084275808340593191)*b0; |
---|
| 1732 | rkE[2] = (KPP_REAL).2870858974881495709929780405075939*b0; |
---|
| 1733 | rkE[3] = (KPP_REAL)(-.4558086256248162565397206448274867e-01)*b0; |
---|
| 1734 | |
---|
| 1735 | /* Sdirk error estimator */ |
---|
| 1736 | rkBgam[0] = ZERO; |
---|
| 1737 | rkBgam[1] = (KPP_REAL).2373339543355109188382583162660537; |
---|
| 1738 | rkBgam[2] = (KPP_REAL).5879873931885192299409334646982414; |
---|
| 1739 | rkBgam[3] = (KPP_REAL)(-.4063577064014232702392531134499046e-01); |
---|
| 1740 | rkBgam[4] = (KPP_REAL).2153144231161121782447335303806955; |
---|
| 1741 | |
---|
| 1742 | /* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */ |
---|
| 1743 | rkTheta[0] = (KPP_REAL)(-2.594040933093095272574031876464493); |
---|
| 1744 | rkTheta[1] = (KPP_REAL)1.824611539036311947589425112250199; |
---|
| 1745 | rkTheta[2] = (KPP_REAL).1856563166634371860478043996459493; |
---|
| 1746 | |
---|
| 1747 | /* ELO = local order of classical error estimator */ |
---|
| 1748 | rkELO = (KPP_REAL)4.0; |
---|
| 1749 | |
---|
| 1750 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1751 | ~~~> Diagonalize the RK matrix: |
---|
| 1752 | rkTinv * inv(rkA) * rkT = |
---|
| 1753 | | rkGamma 0 0 | |
---|
| 1754 | | 0 rkAlpha -rkBeta | |
---|
| 1755 | | 0 rkBeta rkAlpha | |
---|
| 1756 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1757 | |
---|
| 1758 | rkGamma = (KPP_REAL)4.644370709252171185822941421408064; |
---|
| 1759 | rkAlpha = (KPP_REAL)3.677814645373914407088529289295970; |
---|
| 1760 | rkBeta = (KPP_REAL)3.508761919567443321903661209182446; |
---|
| 1761 | |
---|
| 1762 | rkT[0][0] = (KPP_REAL)(.7215185205520017032081769924397664e-01); |
---|
| 1763 | rkT[0][1] = (KPP_REAL)(-.8224123057363067064866206597516454e-01); |
---|
| 1764 | rkT[0][2] = (KPP_REAL)(-.6012073861930850173085948921439054e-01); |
---|
| 1765 | rkT[1][0] = (KPP_REAL).1188325787412778070708888193730294; |
---|
| 1766 | rkT[1][1] = (KPP_REAL)(.5306509074206139504614411373957448e-01); |
---|
| 1767 | rkT[1][2] = (KPP_REAL).3162050511322915732224862926182701; |
---|
| 1768 | rkT[2][0] = ONE; |
---|
| 1769 | rkT[2][1] = ONE; |
---|
| 1770 | rkT[2][2] = ZERO; |
---|
| 1771 | |
---|
| 1772 | rkTinv[0][0] = (KPP_REAL)5.991698084937800775649580743981285; |
---|
| 1773 | rkTinv[0][1] = (KPP_REAL)1.139214295155735444567002236934009; |
---|
| 1774 | rkTinv[0][2] = (KPP_REAL).4323121137838583855696375901180497; |
---|
| 1775 | rkTinv[1][0] = (KPP_REAL)(-5.991698084937800775649580743981285); |
---|
| 1776 | rkTinv[1][1] = (KPP_REAL)(-1.139214295155735444567002236934009); |
---|
| 1777 | rkTinv[1][2] = (KPP_REAL).5676878862161416144303624098819503; |
---|
| 1778 | rkTinv[2][0] = (KPP_REAL)(-1.246213273586231410815571640493082); |
---|
| 1779 | rkTinv[2][1] = (KPP_REAL)2.925559646192313662599230367054972; |
---|
| 1780 | rkTinv[2][2] = (KPP_REAL)(-.2577352012734324923468722836888244); |
---|
| 1781 | |
---|
| 1782 | rkTinvAinv[0][0] = (KPP_REAL)27.82766708436744962047620566703329; |
---|
| 1783 | rkTinvAinv[0][1] = (KPP_REAL)5.290933503982655311815946575100597; |
---|
| 1784 | rkTinvAinv[0][2] = (KPP_REAL)2.007817718512643701322151051660114; |
---|
| 1785 | rkTinvAinv[1][0] = (KPP_REAL)(-17.66368928942422710690385180065675); |
---|
| 1786 | rkTinvAinv[1][1] = (KPP_REAL)(-14.45491129892587782538830044147713); |
---|
| 1787 | rkTinvAinv[1][2] = (KPP_REAL)2.992182281487356298677848948339886; |
---|
| 1788 | rkTinvAinv[2][0] = (KPP_REAL)(-25.60678350282974256072419392007303); |
---|
| 1789 | rkTinvAinv[2][1] = (KPP_REAL)6.762434375611708328910623303779923; |
---|
| 1790 | rkTinvAinv[2][2] = (KPP_REAL)1.043979339483109825041215970036771; |
---|
| 1791 | |
---|
| 1792 | rkAinvT[0][0] = (KPP_REAL).3350999483034677402618981153470483; |
---|
| 1793 | rkAinvT[0][1] = (KPP_REAL)(-.5134173605009692329246186488441294); |
---|
| 1794 | rkAinvT[0][2] = (KPP_REAL)(.6745196507033116204327635673208923e-01); |
---|
| 1795 | rkAinvT[1][0] = (KPP_REAL).5519025480108928886873752035738885; |
---|
| 1796 | rkAinvT[1][1] = (KPP_REAL)1.304651810077110066076640761092008; |
---|
| 1797 | rkAinvT[1][2] = (KPP_REAL).9767507983414134987545585703726984; |
---|
| 1798 | rkAinvT[2][0] = (KPP_REAL)4.644370709252171185822941421408064; |
---|
| 1799 | rkAinvT[2][1] = (KPP_REAL)3.677814645373914407088529289295970; |
---|
| 1800 | rkAinvT[2][2] = (KPP_REAL)(-3.508761919567443321903661209182446); |
---|
| 1801 | } /* Gauss_Coefficients */ |
---|
| 1802 | |
---|
| 1803 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1804 | The coefficients of the 3-stage Gauss method |
---|
| 1805 | (given to ~30 accurate digits) |
---|
| 1806 | The parameter b3 can be chosen by the user |
---|
| 1807 | to tune the error estimator |
---|
| 1808 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1809 | void Radau1A_Coefficients() |
---|
| 1810 | { |
---|
| 1811 | /* The coefficients of the Radau1A method */ |
---|
| 1812 | KPP_REAL b0; |
---|
| 1813 | /*b0 = (KPP_REAL)0.3; */ |
---|
| 1814 | b0 = (KPP_REAL)0.1; |
---|
| 1815 | |
---|
| 1816 | /* The coefficients of the Radau1A method */ |
---|
| 1817 | rkMethod = R1A; |
---|
| 1818 | |
---|
| 1819 | rkA[0][0] = (KPP_REAL) .1111111111111111111111111111111111; |
---|
| 1820 | rkA[0][1] = (KPP_REAL)(-.1916383190435098943442935597058829); |
---|
| 1821 | rkA[0][2] = (KPP_REAL)(.8052720793239878323318244859477174e-01); |
---|
| 1822 | rkA[1][0] = (KPP_REAL).1111111111111111111111111111111111; |
---|
| 1823 | rkA[1][1] = (KPP_REAL).2920734116652284630205027458970589; |
---|
| 1824 | rkA[1][2] = (KPP_REAL)(-.481334970546573839513422644787591e-01); |
---|
| 1825 | rkA[2][0] = (KPP_REAL).1111111111111111111111111111111111; |
---|
| 1826 | rkA[2][1] = (KPP_REAL).5370223859435462728402311533676479; |
---|
| 1827 | rkA[2][2] = (KPP_REAL).1968154772236604258683861429918299; |
---|
| 1828 | |
---|
| 1829 | rkB[0] = (KPP_REAL).1111111111111111111111111111111111; |
---|
| 1830 | rkB[1] = (KPP_REAL).5124858261884216138388134465196080; |
---|
| 1831 | rkB[2] = (KPP_REAL).3764030627004672750500754423692808; |
---|
| 1832 | |
---|
| 1833 | rkC[0] = ZERO; |
---|
| 1834 | rkC[1] = (KPP_REAL).3550510257216821901802715925294109; |
---|
| 1835 | rkC[2] = (KPP_REAL).8449489742783178098197284074705891; |
---|
| 1836 | |
---|
| 1837 | /* Classical error estimator, embedded solution: */ |
---|
| 1838 | rkBhat[0] = b0; |
---|
| 1839 | rkBhat[1] = (KPP_REAL).11111111111111111111111111111111111-b0; |
---|
| 1840 | rkBhat[2] = (KPP_REAL).51248582618842161383881344651960810; |
---|
| 1841 | rkBhat[3] = (KPP_REAL).37640306270046727505007544236928079; |
---|
| 1842 | |
---|
| 1843 | /* New solution: H* Sum B_j*f(Z_j) = Sum D_j*Z_j */ |
---|
| 1844 | rkD[0] = (KPP_REAL).3333333333333333333333333333333333; |
---|
| 1845 | rkD[1] = (KPP_REAL)(-.8914115380582557157653087040196127); |
---|
| 1846 | rkD[2] = (KPP_REAL)(1.558078204724922382431975370686279); |
---|
| 1847 | |
---|
| 1848 | /* Classical error estimator: */ |
---|
| 1849 | /* H* Sum (b_j-bhat_j) f(Z_j) = H*E(0)*F(0) + Sum E_j Z_j */ |
---|
| 1850 | rkE[0] = (KPP_REAL).2748888295956773677478286035994148*b0; |
---|
| 1851 | rkE[1] = (KPP_REAL)(-1.374444147978386838739143017997074)*b0; |
---|
| 1852 | rkE[2] = (KPP_REAL)(-1.335337922441686804550326197041126)*b0; |
---|
| 1853 | rkE[3] = (KPP_REAL).235782604058977333559011782643466*b0; |
---|
| 1854 | |
---|
| 1855 | /* Sdirk error estimator */ |
---|
| 1856 | rkBgam[0] = ZERO; |
---|
| 1857 | rkBgam[1] = (KPP_REAL)(.1948150124588532186183490991130616e-01); |
---|
| 1858 | rkBgam[2] = (KPP_REAL).7575249005733381398986810981093584; |
---|
| 1859 | rkBgam[3] = (KPP_REAL)(-.518952314149008295083446116200793e-01); |
---|
| 1860 | rkBgam[4] = (KPP_REAL).2748888295956773677478286035994148; |
---|
| 1861 | |
---|
| 1862 | /* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */ |
---|
| 1863 | rkTheta[0] = (KPP_REAL)(-1.224370034375505083904362087063351); |
---|
| 1864 | rkTheta[1] = (KPP_REAL).9340045331532641409047527962010133; |
---|
| 1865 | rkTheta[2] = (KPP_REAL).4656990124352088397561234800640929; |
---|
| 1866 | |
---|
| 1867 | /* ELO = local order of classical error estimator */ |
---|
| 1868 | rkELO = (KPP_REAL)4.0; |
---|
| 1869 | |
---|
| 1870 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1871 | ~~~> Diagonalize the RK matrix: |
---|
| 1872 | rkTinv * inv(rkA) * rkT = |
---|
| 1873 | | rkGamma 0 0 | |
---|
| 1874 | | 0 rkAlpha -rkBeta | |
---|
| 1875 | | 0 rkBeta rkAlpha | |
---|
| 1876 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1877 | |
---|
| 1878 | rkGamma = (KPP_REAL)3.637834252744495732208418513577775; |
---|
| 1879 | rkAlpha = (KPP_REAL)2.681082873627752133895790743211112; |
---|
| 1880 | rkBeta = (KPP_REAL)3.050430199247410569426377624787569; |
---|
| 1881 | |
---|
| 1882 | rkT[0][0] = (KPP_REAL).424293819848497965354371036408369; |
---|
| 1883 | rkT[0][1] = (KPP_REAL)(-.3235571519651980681202894497035503); |
---|
| 1884 | rkT[0][2] = (KPP_REAL)(-.522137786846287839586599927945048); |
---|
| 1885 | rkT[1][0] = (KPP_REAL)(.57594609499806128896291585429339e-01); |
---|
| 1886 | rkT[1][1] = (KPP_REAL)(.3148663231849760131614374283783e-02); |
---|
| 1887 | rkT[1][2] = (KPP_REAL).452429247674359778577728510381731; |
---|
| 1888 | rkT[2][0] = ONE; |
---|
| 1889 | rkT[2][1] = ONE; |
---|
| 1890 | rkT[2][2] = ZERO; |
---|
| 1891 | |
---|
| 1892 | rkTinv[0][0] = (KPP_REAL)1.233523612685027760114769983066164; |
---|
| 1893 | rkTinv[0][1] = (KPP_REAL)1.423580134265707095505388133369554; |
---|
| 1894 | rkTinv[0][2] = (KPP_REAL).3946330125758354736049045150429624; |
---|
| 1895 | rkTinv[1][0] = (KPP_REAL)(-1.233523612685027760114769983066164); |
---|
| 1896 | rkTinv[1][1] = (KPP_REAL)(-1.423580134265707095505388133369554); |
---|
| 1897 | rkTinv[1][2] = (KPP_REAL).6053669874241645263950954849570376; |
---|
| 1898 | rkTinv[2][0] = (KPP_REAL)(-.1484438963257383124456490049673414); |
---|
| 1899 | rkTinv[2][1] = (KPP_REAL)2.038974794939896109682070471785315; |
---|
| 1900 | rkTinv[2][2] = (KPP_REAL)(-.544501292892686735299355831692542e-01); |
---|
| 1901 | |
---|
| 1902 | rkTinvAinv[0][0] = (KPP_REAL)4.487354449794728738538663081025420; |
---|
| 1903 | rkTinvAinv[0][1] = (KPP_REAL)5.178748573958397475446442544234494; |
---|
| 1904 | rkTinvAinv[0][2] = (KPP_REAL)1.435609490412123627047824222335563; |
---|
| 1905 | rkTinvAinv[1][0] = (KPP_REAL)(-2.854361287939276673073807031221493); |
---|
| 1906 | rkTinvAinv[1][1] = (KPP_REAL)(-1.003648660720543859000994063139137e+01); |
---|
| 1907 | rkTinvAinv[1][2] = (KPP_REAL)1.789135380979465422050817815017383; |
---|
| 1908 | rkTinvAinv[2][0] = (KPP_REAL)(-4.160768067752685525282947313530352); |
---|
| 1909 | rkTinvAinv[2][1] = (KPP_REAL)1.124128569859216916690209918405860; |
---|
| 1910 | rkTinvAinv[2][2] = (KPP_REAL)1.700644430961823796581896350418417; |
---|
| 1911 | |
---|
| 1912 | rkAinvT[0][0] = (KPP_REAL)1.543510591072668287198054583233180; |
---|
| 1913 | rkAinvT[0][1] = (KPP_REAL)(-2.460228411937788329157493833295004); |
---|
| 1914 | rkAinvT[0][2] = (KPP_REAL)(-.412906170450356277003910443520499); |
---|
| 1915 | rkAinvT[1][0] = (KPP_REAL).209519643211838264029272585946993; |
---|
| 1916 | rkAinvT[1][1] = (KPP_REAL)1.388545667194387164417459732995766; |
---|
| 1917 | rkAinvT[1][2] = (KPP_REAL)1.20339553005832004974976023130002; |
---|
| 1918 | rkAinvT[2][0] = (KPP_REAL)3.637834252744495732208418513577775; |
---|
| 1919 | rkAinvT[2][1] = (KPP_REAL)2.681082873627752133895790743211112; |
---|
| 1920 | rkAinvT[2][2] = (KPP_REAL)(-3.050430199247410569426377624787569); |
---|
| 1921 | } /* Radau1A_Coefficients */ |
---|
| 1922 | |
---|
| 1923 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1924 | The coefficients of the 4-stage Lobatto-3A method |
---|
| 1925 | (given to ~30 accurate digits) |
---|
| 1926 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1927 | void Lobatto3A_Coefficients() |
---|
| 1928 | { |
---|
| 1929 | /* The coefficients of the Lobatto-3A method */ |
---|
| 1930 | rkMethod = L3A; |
---|
| 1931 | |
---|
| 1932 | rkA[0][0] = ZERO; |
---|
| 1933 | rkA[0][1] = ZERO; |
---|
| 1934 | rkA[0][2] = ZERO; |
---|
| 1935 | rkA[0][3] = ZERO; |
---|
| 1936 | rkA[1][0] = (KPP_REAL).11030056647916491413674311390609397; |
---|
| 1937 | rkA[1][1] = (KPP_REAL).1896994335208350858632568860939060; |
---|
| 1938 | rkA[1][2] = (KPP_REAL)(-.339073642291438837776604807792215e-01); |
---|
| 1939 | rkA[1][3] = (KPP_REAL)(.1030056647916491413674311390609397e-01); |
---|
| 1940 | rkA[2][0] = (KPP_REAL)(.73032766854168419196590219427239365e-01); |
---|
| 1941 | rkA[2][1] = (KPP_REAL).4505740308958105504443271474458881; |
---|
| 1942 | rkA[2][2] = (KPP_REAL).2269672331458315808034097805727606; |
---|
| 1943 | rkA[2][3] = (KPP_REAL)(-.2696723314583158080340978057276063e-01); |
---|
| 1944 | rkA[3][0] = (KPP_REAL)(.83333333333333333333333333333333333e-01); |
---|
| 1945 | rkA[3][1] = (KPP_REAL).4166666666666666666666666666666667; |
---|
| 1946 | rkA[3][2] = (KPP_REAL).4166666666666666666666666666666667; |
---|
| 1947 | rkA[3][3] = (KPP_REAL)(.8333333333333333333333333333333333e-01); |
---|
| 1948 | |
---|
| 1949 | rkB[0] = (KPP_REAL)(.83333333333333333333333333333333333e-01); |
---|
| 1950 | rkB[1] = (KPP_REAL).4166666666666666666666666666666667; |
---|
| 1951 | rkB[2] = (KPP_REAL).4166666666666666666666666666666667; |
---|
| 1952 | rkB[3] = (KPP_REAL)(.8333333333333333333333333333333333e-01); |
---|
| 1953 | |
---|
| 1954 | rkC[0] = ZERO; |
---|
| 1955 | rkC[1] = (KPP_REAL).2763932022500210303590826331268724; |
---|
| 1956 | rkC[2] = (KPP_REAL).7236067977499789696409173668731276; |
---|
| 1957 | rkC[3] = ONE; |
---|
| 1958 | |
---|
| 1959 | /* New solution: H*Sum B_j*f(Z_j) = Sum D_j*Z_j */ |
---|
| 1960 | rkD[0] = ZERO; |
---|
| 1961 | rkD[1] = ZERO; |
---|
| 1962 | rkD[2] = ZERO; |
---|
| 1963 | rkD[3] = ONE; |
---|
| 1964 | |
---|
| 1965 | /* Classical error estimator, embedded solution: */ |
---|
| 1966 | rkBhat[0] = (KPP_REAL)(.90909090909090909090909090909090909e-01); |
---|
| 1967 | rkBhat[1] = (KPP_REAL).39972675774621371442114262372173276; |
---|
| 1968 | rkBhat[2] = (KPP_REAL).43360657558711961891219070961160058; |
---|
| 1969 | rkBhat[3] = (KPP_REAL)(.15151515151515151515151515151515152e-01); |
---|
| 1970 | |
---|
| 1971 | /* Classical error estimator: */ |
---|
| 1972 | /* H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j */ |
---|
| 1973 | |
---|
| 1974 | rkE[0] = (KPP_REAL)(.1957403846510110711315759367097231e-01); |
---|
| 1975 | rkE[1] = (KPP_REAL)(-.1986820345632580910316020806676438); |
---|
| 1976 | rkE[2] = (KPP_REAL).1660586371214229125096727578826900; |
---|
| 1977 | rkE[3] = (KPP_REAL)(-.9787019232550553556578796835486154e-01); |
---|
| 1978 | |
---|
| 1979 | /* Sdirk error estimator: */ |
---|
| 1980 | rkF[0] = ZERO; |
---|
| 1981 | rkF[1] = (KPP_REAL)(-.66535815876916686607437314126436349); |
---|
| 1982 | rkF[2] = (KPP_REAL)1.7419302743497277572980407931678409; |
---|
| 1983 | rkF[3] = (KPP_REAL)(-1.2918865386966730694684011822841728); |
---|
| 1984 | |
---|
| 1985 | /* ELO = local order of classical error estimator */ |
---|
| 1986 | rkELO = (KPP_REAL)4.0; |
---|
| 1987 | |
---|
| 1988 | /* Sdirk error estimator: */ |
---|
| 1989 | rkBgam[0] = (KPP_REAL)(.2950472755430528877214995073815946e-01); |
---|
| 1990 | rkBgam[1] = (KPP_REAL).5370310883226113978352873633882769; |
---|
| 1991 | rkBgam[2] = (KPP_REAL).2963022450107219354980459699450564; |
---|
| 1992 | rkBgam[3] = (KPP_REAL)(-.7815248400375080035021681445218837e-01); |
---|
| 1993 | rkBgam[4] = (KPP_REAL).2153144231161121782447335303806956; |
---|
| 1994 | |
---|
| 1995 | /* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */ |
---|
| 1996 | rkTheta[0] = (KPP_REAL)0.0; |
---|
| 1997 | rkTheta[1] = (KPP_REAL)(-.6653581587691668660743731412643631); |
---|
| 1998 | rkTheta[2] = (KPP_REAL)1.741930274349727757298040793167842; |
---|
| 1999 | rkTheta[3] = (KPP_REAL)(-.291886538696673069468401182284174); |
---|
| 2000 | |
---|
| 2001 | |
---|
| 2002 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 2003 | ~~~> Diagonalize the RK matrix: |
---|
| 2004 | rkTinv * inv(rkA) * rkT = |
---|
| 2005 | | rkGamma 0 0 | |
---|
| 2006 | | 0 rkAlpha -rkBeta | |
---|
| 2007 | | 0 rkBeta rkAlpha | |
---|
| 2008 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2009 | |
---|
| 2010 | rkGamma = (KPP_REAL)4.644370709252171185822941421408063; |
---|
| 2011 | rkAlpha = (KPP_REAL)3.677814645373914407088529289295968; |
---|
| 2012 | rkBeta = (KPP_REAL)3.508761919567443321903661209182446; |
---|
| 2013 | |
---|
| 2014 | rkT[0][0] = (KPP_REAL)(.5303036326129938105898786144870856e-01); |
---|
| 2015 | rkT[0][1] = (KPP_REAL)(-.7776129960563076320631956091016914e-01); |
---|
| 2016 | rkT[0][2] = (KPP_REAL)(.6043307469475508514468017399717112e-02); |
---|
| 2017 | rkT[1][0] = (KPP_REAL).2637242522173698467283726114649606; |
---|
| 2018 | rkT[1][1] = (KPP_REAL).2193839918662961493126393244533346; |
---|
| 2019 | rkT[1][2] = (KPP_REAL).3198765142300936188514264752235344; |
---|
| 2020 | rkT[2][0] = ONE; |
---|
| 2021 | rkT[2][1] = ZERO; |
---|
| 2022 | rkT[2][2] = ZERO; |
---|
| 2023 | |
---|
| 2024 | rkTinv[0][0] = (KPP_REAL)7.695032983257654470769069079238553; |
---|
| 2025 | rkTinv[0][1] = (KPP_REAL)(-.1453793830957233720334601186354032); |
---|
| 2026 | rkTinv[0][2] = (KPP_REAL).6302696746849084900422461036874826; |
---|
| 2027 | rkTinv[1][0] = (KPP_REAL)(-7.695032983257654470769069079238553); |
---|
| 2028 | rkTinv[1][1] = (KPP_REAL).1453793830957233720334601186354032; |
---|
| 2029 | rkTinv[1][2] = (KPP_REAL).3697303253150915099577538963125174; |
---|
| 2030 | rkTinv[2][0] = (KPP_REAL)(-1.066660885401270392058552736086173); |
---|
| 2031 | rkTinv[2][1] = (KPP_REAL)3.146358406832537460764521760668932; |
---|
| 2032 | rkTinv[2][2] = (KPP_REAL)(-.7732056038202974770406168510664738); |
---|
| 2033 | |
---|
| 2034 | rkTinvAinv[0][0] = (KPP_REAL)35.73858579417120341641749040405149; |
---|
| 2035 | rkTinvAinv[0][1] = (KPP_REAL)(-.675195748578927863668368190236025); |
---|
| 2036 | rkTinvAinv[0][2] = (KPP_REAL)2.927206016036483646751158874041632; |
---|
| 2037 | rkTinvAinv[1][0] = (KPP_REAL)(-24.55824590667225493437162206039511); |
---|
| 2038 | rkTinvAinv[1][1] = (KPP_REAL)(-10.50514413892002061837750015342036); |
---|
| 2039 | rkTinvAinv[1][2] = (KPP_REAL)4.072793983963516353248841125958369; |
---|
| 2040 | rkTinvAinv[2][0] = (KPP_REAL)(-30.92301972744621647251975054630589); |
---|
| 2041 | rkTinvAinv[2][1] = (KPP_REAL)12.08182467154052413351908559269928; |
---|
| 2042 | rkTinvAinv[2][2] = (KPP_REAL)(-1.546411207640594954081233702132946); |
---|
| 2043 | |
---|
| 2044 | rkAinvT[0][0] = (KPP_REAL).2462926658317812882584158369803835; |
---|
| 2045 | rkAinvT[0][1] = (KPP_REAL)(-.2647871194157644619747121197289574); |
---|
| 2046 | rkAinvT[0][2] = (KPP_REAL).2950720515900466654896406799284586; |
---|
| 2047 | rkAinvT[1][0] = (KPP_REAL)1.224833192317784474576995878738004; |
---|
| 2048 | rkAinvT[1][1] = (KPP_REAL)1.929224190340981580557006261869763; |
---|
| 2049 | rkAinvT[1][2] = (KPP_REAL).4066803323234419988910915619080306; |
---|
| 2050 | rkAinvT[2][0] = (KPP_REAL)4.644370709252171185822941421408064; |
---|
| 2051 | rkAinvT[2][1] = (KPP_REAL)3.677814645373914407088529289295968; |
---|
| 2052 | rkAinvT[2][2] = (KPP_REAL)(-3.508761919567443321903661209182446); |
---|
| 2053 | } /* Lobatto3A_Coefficients */ |
---|
| 2054 | |
---|
| 2055 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 2056 | END SUBROUTINE RungeKutta ! and all its internal procedures |
---|
| 2057 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2058 | |
---|
| 2059 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2060 | void FUN_CHEM(KPP_REAL T, KPP_REAL V[], KPP_REAL FCT[]) |
---|
| 2061 | { |
---|
| 2062 | KPP_REAL Told; |
---|
| 2063 | |
---|
| 2064 | Told = TIME; |
---|
| 2065 | TIME = T; |
---|
| 2066 | Update_SUN(); |
---|
| 2067 | Update_RCONST(); |
---|
| 2068 | Update_PHOTO(); |
---|
| 2069 | TIME = Told; |
---|
| 2070 | |
---|
| 2071 | Fun(V, FIX, RCONST, FCT); |
---|
| 2072 | |
---|
| 2073 | } /* FUN_CHEM */ |
---|
| 2074 | |
---|
| 2075 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 2076 | void JAC_CHEM (KPP_REAL T, KPP_REAL V[], KPP_REAL JF[]) |
---|
| 2077 | { |
---|
| 2078 | KPP_REAL Told; |
---|
| 2079 | |
---|
| 2080 | Told = TIME; |
---|
| 2081 | TIME = T; |
---|
| 2082 | Update_SUN(); |
---|
| 2083 | Update_RCONST(); |
---|
| 2084 | Update_PHOTO(); |
---|
| 2085 | TIME = Told; |
---|
| 2086 | Jac_SP(V, FIX, RCONST, JF); |
---|
| 2087 | } /* JAC_CHEM */ |
---|