[2696] | 1 | |
---|
| 2 | #define MAX(a,b) ( ((a) >= (b)) ?(a):(b) ) |
---|
| 3 | #define MIN(b,c) ( ((b) < (c)) ?(b):(c) ) |
---|
| 4 | #define ABS(x) ( ((x) >= 0 ) ?(x):(-x) ) |
---|
| 5 | #define SQRT(d) ( pow((d),0.5) ) |
---|
| 6 | #define SIGN(x) ( ((x) >= 0 ) ?[0]:(-1) ) |
---|
| 7 | |
---|
| 8 | /*~~> Numerical constants */ |
---|
| 9 | #define ZERO (KPP_REAL)0.0 |
---|
| 10 | #define ONE (KPP_REAL)1.0 |
---|
| 11 | #define HALF (KPP_REAL)0.5 |
---|
| 12 | #define DeltaMin (KPP_REAL)1.0e-6 |
---|
| 13 | |
---|
| 14 | /*~~~> Collect statistics: global variables */ |
---|
| 15 | int Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng; |
---|
| 16 | |
---|
| 17 | |
---|
| 18 | /*~~~> Function headers */ |
---|
| 19 | void FunTemplate(KPP_REAL, KPP_REAL [], KPP_REAL []); |
---|
| 20 | void JacTemplate(KPP_REAL, KPP_REAL [], KPP_REAL []) ; |
---|
| 21 | int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, |
---|
| 22 | KPP_REAL AbsTol[], KPP_REAL RelTol[], |
---|
| 23 | void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), |
---|
| 24 | void (*ode_Jac)(KPP_REAL, KPP_REAL [], KPP_REAL []), |
---|
| 25 | KPP_REAL RPAR[], int IPAR[]); |
---|
| 26 | int RosenbrockIntegrator( |
---|
| 27 | KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend , |
---|
| 28 | KPP_REAL AbsTol[], KPP_REAL RelTol[], |
---|
| 29 | void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), |
---|
| 30 | void (*ode_Jac)(KPP_REAL, KPP_REAL [], KPP_REAL []), |
---|
| 31 | int ros_S, |
---|
| 32 | KPP_REAL ros_M[], KPP_REAL ros_E[], |
---|
| 33 | KPP_REAL ros_A[], KPP_REAL ros_C[], |
---|
| 34 | KPP_REAL ros_Alpha[],KPP_REAL ros_Gamma[], |
---|
| 35 | KPP_REAL ros_ELO, char ros_NewF[], |
---|
| 36 | char Autonomous, char VectorTol, int Max_no_steps, |
---|
| 37 | KPP_REAL Roundoff, KPP_REAL Hmin, KPP_REAL Hmax, KPP_REAL Hstart, |
---|
| 38 | KPP_REAL FacMin, KPP_REAL FacMax, KPP_REAL FacRej, KPP_REAL FacSafe, |
---|
| 39 | KPP_REAL *Texit, KPP_REAL *Hexit ); |
---|
| 40 | char ros_PrepareMatrix ( |
---|
| 41 | KPP_REAL* H, |
---|
| 42 | int Direction, KPP_REAL gam, KPP_REAL Jac0[], |
---|
| 43 | KPP_REAL Ghimj[], int Pivot[] ); |
---|
| 44 | KPP_REAL ros_ErrorNorm ( |
---|
| 45 | KPP_REAL Y[], KPP_REAL Ynew[], KPP_REAL Yerr[], |
---|
| 46 | KPP_REAL AbsTol[], KPP_REAL RelTol[], |
---|
| 47 | char VectorTol ); |
---|
| 48 | int ros_ErrorMsg(int Code, KPP_REAL T, KPP_REAL H); |
---|
| 49 | void ros_FunTimeDerivative ( |
---|
| 50 | KPP_REAL T, KPP_REAL Roundoff, |
---|
| 51 | KPP_REAL Y[], KPP_REAL Fcn0[], |
---|
| 52 | void ode_Fun(KPP_REAL, KPP_REAL [], KPP_REAL []), |
---|
| 53 | KPP_REAL dFdT[] ); |
---|
| 54 | void Fun( KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[] ); |
---|
| 55 | void Jac_SP( KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[] ); |
---|
| 56 | void FunTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Ydot[] ); |
---|
| 57 | void JacTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Ydot[] ); |
---|
| 58 | void DecompTemplate( KPP_REAL A[], int Pivot[], int* ising ); |
---|
| 59 | void SolveTemplate( KPP_REAL A[], int Pivot[], KPP_REAL b[] ); |
---|
| 60 | void WCOPY(int N, KPP_REAL X[], int incX, KPP_REAL Y[], int incY); |
---|
| 61 | void WAXPY(int N, KPP_REAL Alpha, KPP_REAL X[], int incX, KPP_REAL Y[], int incY ); |
---|
| 62 | void WSCAL(int N, KPP_REAL Alpha, KPP_REAL X[], int incX); |
---|
| 63 | KPP_REAL WLAMCH( char C ); |
---|
| 64 | void Ros2 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], |
---|
| 65 | KPP_REAL ros_M[], KPP_REAL ros_E[], |
---|
| 66 | KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], |
---|
| 67 | char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ); |
---|
| 68 | void Ros3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], |
---|
| 69 | KPP_REAL ros_M[], KPP_REAL ros_E[], |
---|
| 70 | KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], |
---|
| 71 | char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ); |
---|
| 72 | void Ros4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], |
---|
| 73 | KPP_REAL ros_M[], KPP_REAL ros_E[], |
---|
| 74 | KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], |
---|
| 75 | char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ); |
---|
| 76 | void Rodas3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], |
---|
| 77 | KPP_REAL ros_M[], KPP_REAL ros_E[], |
---|
| 78 | KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], |
---|
| 79 | char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ); |
---|
| 80 | void Rodas4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], |
---|
| 81 | KPP_REAL ros_M[], KPP_REAL ros_E[], |
---|
| 82 | KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], |
---|
| 83 | char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ); |
---|
| 84 | int KppDecomp( KPP_REAL A[] ); |
---|
| 85 | void KppSolve ( KPP_REAL A[], KPP_REAL b[] ); |
---|
| 86 | void Update_SUN(); |
---|
| 87 | void Update_RCONST(); |
---|
| 88 | |
---|
| 89 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 90 | void INTEGRATE( KPP_REAL TIN, KPP_REAL TOUT ) |
---|
| 91 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 92 | { |
---|
| 93 | static KPP_REAL RPAR[20]; |
---|
| 94 | static int i, IERR, IPAR[20]; |
---|
| 95 | static int Ns=0, Na=0, Nr=0, Ng=0; |
---|
| 96 | |
---|
| 97 | for ( i = 0; i < 20; i++ ) { |
---|
| 98 | IPAR[i] = 0; |
---|
| 99 | RPAR[i] = ZERO; |
---|
| 100 | } /* for */ |
---|
| 101 | |
---|
| 102 | |
---|
| 103 | IPAR[0] = 0; /* non-autonomous */ |
---|
| 104 | IPAR[1] = 1; /* vector tolerances */ |
---|
| 105 | RPAR[2] = STEPMIN; /* starting step */ |
---|
| 106 | IPAR[3] = 5; /* choice of the method */ |
---|
| 107 | |
---|
| 108 | IERR = Rosenbrock(VAR, TIN, TOUT, |
---|
| 109 | ATOL, RTOL, |
---|
| 110 | &FunTemplate, &JacTemplate, |
---|
| 111 | RPAR, IPAR); |
---|
| 112 | |
---|
| 113 | |
---|
| 114 | Ns=Ns+IPAR[12]; |
---|
| 115 | Na=Na+IPAR[13]; |
---|
| 116 | Nr=Nr+IPAR[14]; |
---|
| 117 | Ng=Ng+IPAR[17]; |
---|
| 118 | printf("\n Step=%d Acc=%d Rej=%d Singular=%d\n", |
---|
| 119 | Ns,Na,Nr,Ng); |
---|
| 120 | |
---|
| 121 | |
---|
| 122 | if (IERR < 0) |
---|
| 123 | printf("\n Rosenbrock: Unsucessful step at T=%g: IERR=%d\n", |
---|
| 124 | TIN,IERR); |
---|
| 125 | |
---|
| 126 | TIN = RPAR[10]; /* Exit time */ |
---|
| 127 | STEPMIN = RPAR[11]; /* Last step */ |
---|
| 128 | |
---|
| 129 | } /* INTEGRATE */ |
---|
| 130 | |
---|
| 131 | |
---|
| 132 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 133 | int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, |
---|
| 134 | KPP_REAL AbsTol[], KPP_REAL RelTol[], |
---|
| 135 | void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), |
---|
| 136 | void (*ode_Jac)(KPP_REAL, KPP_REAL [], KPP_REAL []), |
---|
| 137 | KPP_REAL RPAR[], int IPAR[]) |
---|
| 138 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 139 | |
---|
| 140 | Solves the system y'=F(t,y) using a Rosenbrock method defined by: |
---|
| 141 | |
---|
| 142 | G = 1/(H*gamma[0]) - ode_Jac(t0,Y0) |
---|
| 143 | T_i = t0 + Alpha(i)*H |
---|
| 144 | Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j |
---|
| 145 | G * K_i = ode_Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j + |
---|
| 146 | gamma(i)*dF/dT(t0, Y0) |
---|
| 147 | Y1 = Y0 + \sum_{j=1}^S M(j)*K_j |
---|
| 148 | |
---|
| 149 | For details on Rosenbrock methods and their implementation consult: |
---|
| 150 | E. Hairer and G. Wanner |
---|
| 151 | "Solving ODEs II. Stiff and differential-algebraic problems". |
---|
| 152 | Springer series in computational mathematics, Springer-Verlag, 1996. |
---|
| 153 | The codes contained in the book inspired this implementation. |
---|
| 154 | |
---|
| 155 | (C) Adrian Sandu, August 2004 |
---|
| 156 | Virginia Polytechnic Institute and State University |
---|
| 157 | Contact: sandu@cs.vt.edu |
---|
| 158 | This implementation is part of KPP - the Kinetic PreProcessor |
---|
| 159 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 160 | |
---|
| 161 | *~~~> INPUT ARGUMENTS: |
---|
| 162 | |
---|
| 163 | - Y(NVAR) = vector of initial conditions (at T=Tstart) |
---|
| 164 | - [Tstart,Tend] = time range of integration |
---|
| 165 | (if Tstart>Tend the integration is performed backwards in time) |
---|
| 166 | - RelTol, AbsTol = user precribed accuracy |
---|
| 167 | - void ode_Fun( T, Y, Ydot ) = ODE function, |
---|
| 168 | returns Ydot = Y' = F(T,Y) |
---|
| 169 | - void ode_Fun( T, Y, Ydot ) = Jacobian of the ODE function, |
---|
| 170 | returns Jcb = dF/dY |
---|
| 171 | - IPAR(1:10) = int inputs parameters |
---|
| 172 | - RPAR(1:10) = real inputs parameters |
---|
| 173 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 174 | |
---|
| 175 | *~~~> OUTPUT ARGUMENTS: |
---|
| 176 | |
---|
| 177 | - Y(NVAR) -> vector of final states (at T->Tend) |
---|
| 178 | - IPAR(11:20) -> int output parameters |
---|
| 179 | - RPAR(11:20) -> real output parameters |
---|
| 180 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 181 | |
---|
| 182 | *~~~> RETURN VALUE (int): |
---|
| 183 | |
---|
| 184 | - IERR -> job status upon return |
---|
| 185 | - succes (positive value) or failure (negative value) - |
---|
| 186 | = 1 : Success |
---|
| 187 | = -1 : Improper value for maximal no of steps |
---|
| 188 | = -2 : Selected Rosenbrock method not implemented |
---|
| 189 | = -3 : Hmin/Hmax/Hstart must be positive |
---|
| 190 | = -4 : FacMin/FacMax/FacRej must be positive |
---|
| 191 | = -5 : Improper tolerance values |
---|
| 192 | = -6 : No of steps exceeds maximum bound |
---|
| 193 | = -7 : Step size too small |
---|
| 194 | = -8 : Matrix is repeatedly singular |
---|
| 195 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 196 | |
---|
| 197 | *~~~> INPUT PARAMETERS: |
---|
| 198 | |
---|
| 199 | Note: For input parameters equal to zero the default values of the |
---|
| 200 | corresponding variables are used. |
---|
| 201 | |
---|
| 202 | IPAR[0] = 1: F = F(y) Independent of T (AUTONOMOUS) |
---|
| 203 | = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS) |
---|
| 204 | IPAR[1] = 0: AbsTol, RelTol are NVAR-dimensional vectors |
---|
| 205 | = 1: AbsTol, RelTol are scalars |
---|
| 206 | IPAR[2] -> maximum number of integration steps |
---|
| 207 | For IPAR[2]=0) the default value of 100000 is used |
---|
| 208 | |
---|
| 209 | IPAR[3] -> selection of a particular Rosenbrock method |
---|
| 210 | = 0 : default method is Rodas3 |
---|
| 211 | = 1 : method is Ros2 |
---|
| 212 | = 2 : method is Ros3 |
---|
| 213 | = 3 : method is Ros4 |
---|
| 214 | = 4 : method is Rodas3 |
---|
| 215 | = 5: method is Rodas4 |
---|
| 216 | |
---|
| 217 | RPAR[0] -> Hmin, lower bound for the integration step size |
---|
| 218 | It is strongly recommended to keep Hmin = ZERO |
---|
| 219 | RPAR[1] -> Hmax, upper bound for the integration step size |
---|
| 220 | RPAR[2] -> Hstart, starting value for the integration step size |
---|
| 221 | |
---|
| 222 | RPAR[3] -> FacMin, lower bound on step decrease factor (default=0.2) |
---|
| 223 | RPAR[4] -> FacMin,upper bound on step increase factor (default=6) |
---|
| 224 | RPAR[5] -> FacRej, step decrease factor after multiple rejections |
---|
| 225 | (default=0.1) |
---|
| 226 | RPAR[6] -> FacSafe, by which the new step is slightly smaller |
---|
| 227 | than the predicted value (default=0.9) |
---|
| 228 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 229 | |
---|
| 230 | *~~~> OUTPUT PARAMETERS: |
---|
| 231 | |
---|
| 232 | Note: each call to Rosenbrock adds the corrent no. of fcn calls |
---|
| 233 | to previous value of IPAR[10], and similar for the other params. |
---|
| 234 | Set IPAR(11:20) = 0 before call to avoid this accumulation. |
---|
| 235 | |
---|
| 236 | IPAR[10] = No. of function calls |
---|
| 237 | IPAR[11] = No. of jacobian calls |
---|
| 238 | IPAR[12] = No. of steps |
---|
| 239 | IPAR[13] = No. of accepted steps |
---|
| 240 | IPAR[14] = No. of rejected steps (except at the beginning) |
---|
| 241 | IPAR[15] = No. of LU decompositions |
---|
| 242 | IPAR[16] = No. of forward/backward substitutions |
---|
| 243 | IPAR[17] = No. of singular matrix decompositions |
---|
| 244 | |
---|
| 245 | RPAR[10] -> Texit, the time corresponding to the |
---|
| 246 | computed Y upon return |
---|
| 247 | RPAR[11] -> Hexit, last accepted step before exit |
---|
| 248 | For multiple restarts, use Hexit as Hstart in the following run |
---|
| 249 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ |
---|
| 250 | { |
---|
| 251 | |
---|
| 252 | /*~~~> The method parameters */ |
---|
| 253 | #define Smax 6 |
---|
| 254 | int Method, ros_S; |
---|
| 255 | KPP_REAL ros_M[Smax], ros_E[Smax]; |
---|
| 256 | KPP_REAL ros_A[Smax*(Smax-1)/2], ros_C[Smax*(Smax-1)/2]; |
---|
| 257 | KPP_REAL ros_Alpha[Smax], ros_Gamma[Smax], ros_ELO; |
---|
| 258 | char ros_NewF[Smax], ros_Name[12]; |
---|
| 259 | /*~~~> Local variables */ |
---|
| 260 | int Max_no_steps, IERR, i, UplimTol; |
---|
| 261 | char Autonomous, VectorTol; |
---|
| 262 | KPP_REAL Roundoff,FacMin,FacMax,FacRej,FacSafe; |
---|
| 263 | KPP_REAL Hmin, Hmax, Hstart, Hexit, Texit; |
---|
| 264 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ |
---|
| 265 | |
---|
| 266 | /*~~~> Initialize statistics */ |
---|
| 267 | Nfun = IPAR[10]; |
---|
| 268 | Njac = IPAR[11]; |
---|
| 269 | Nstp = IPAR[12]; |
---|
| 270 | Nacc = IPAR[13]; |
---|
| 271 | Nrej = IPAR[14]; |
---|
| 272 | Ndec = IPAR[15]; |
---|
| 273 | Nsol = IPAR[16]; |
---|
| 274 | Nsng = IPAR[17]; |
---|
| 275 | |
---|
| 276 | /*~~~> Autonomous or time dependent ODE. Default is time dependent. */ |
---|
| 277 | Autonomous = !(IPAR[0] == 0); |
---|
| 278 | |
---|
| 279 | /*~~~> For Scalar tolerances (IPAR[1] != 0) the code uses AbsTol[0] and RelTol[0] |
---|
| 280 | ! For Vector tolerances (IPAR[1] == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR) */ |
---|
| 281 | if (IPAR[1] == 0) { |
---|
| 282 | VectorTol = 1; UplimTol = KPP_NVAR; |
---|
| 283 | } else { |
---|
| 284 | VectorTol = 0; UplimTol = 1; |
---|
| 285 | } /* end if */ |
---|
| 286 | |
---|
| 287 | /*~~~> The maximum number of steps admitted */ |
---|
| 288 | if (IPAR[2] == 0) |
---|
| 289 | Max_no_steps = 100000; |
---|
| 290 | else |
---|
| 291 | Max_no_steps=IPAR[2]; |
---|
| 292 | if (Max_no_steps < 0) { |
---|
| 293 | printf("\n User-selected max no. of steps: IPAR[2]=%d\n",IPAR[2]); |
---|
| 294 | return ros_ErrorMsg(-1,Tstart,ZERO); |
---|
| 295 | } /* end if */ |
---|
| 296 | |
---|
| 297 | /*~~~> The particular Rosenbrock method chosen */ |
---|
| 298 | if (IPAR[3] == 0) |
---|
| 299 | Method = 3; |
---|
| 300 | else |
---|
| 301 | Method = IPAR[3]; |
---|
| 302 | if ( (IPAR[3] < 1) || (IPAR[3] > 5) ){ |
---|
| 303 | printf("\n User-selected Rosenbrock method: IPAR[3]=%d\n",IPAR[3]); |
---|
| 304 | return ros_ErrorMsg(-2,Tstart,ZERO); |
---|
| 305 | } /* end if */ |
---|
| 306 | |
---|
| 307 | /*~~~> Unit Roundoff (1+Roundoff>1) */ |
---|
| 308 | Roundoff = WLAMCH('E'); |
---|
| 309 | |
---|
| 310 | /*~~~> Lower bound on the step size: (positive value) */ |
---|
| 311 | Hmin = RPAR[0]; |
---|
| 312 | if (RPAR[0] < ZERO) { |
---|
| 313 | printf("\n User-selected Hmin: RPAR[0]=%e\n", RPAR[0]); |
---|
| 314 | return ros_ErrorMsg(-3,Tstart,ZERO); |
---|
| 315 | } /* end if */ |
---|
| 316 | /*~~~> Upper bound on the step size: (positive value) */ |
---|
| 317 | if (RPAR[1] == ZERO) |
---|
| 318 | Hmax = ABS(Tend-Tstart); |
---|
| 319 | else |
---|
| 320 | Hmax = MIN(ABS(RPAR[1]),ABS(Tend-Tstart)); |
---|
| 321 | if (RPAR[1] < ZERO) { |
---|
| 322 | printf("\n User-selected Hmax: RPAR[1]=%e\n", RPAR[1]); |
---|
| 323 | return ros_ErrorMsg(-3,Tstart,ZERO); |
---|
| 324 | } /* end if */ |
---|
| 325 | /*~~~> Starting step size: (positive value) */ |
---|
| 326 | if (RPAR[2] == ZERO) |
---|
| 327 | Hstart = MAX(Hmin,DeltaMin); |
---|
| 328 | else |
---|
| 329 | Hstart = MIN(ABS(RPAR[2]),ABS(Tend-Tstart)); |
---|
| 330 | if (RPAR[2] < ZERO) { |
---|
| 331 | printf("\n User-selected Hstart: RPAR[2]=%e\n", RPAR[2]); |
---|
| 332 | return ros_ErrorMsg(-3,Tstart,ZERO); |
---|
| 333 | } /* end if */ |
---|
| 334 | /*~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax */ |
---|
| 335 | if (RPAR[3] == ZERO) |
---|
| 336 | FacMin = (KPP_REAL)0.2; |
---|
| 337 | else |
---|
| 338 | FacMin = RPAR[3]; |
---|
| 339 | if (RPAR[3] < ZERO) { |
---|
| 340 | printf("\n User-selected FacMin: RPAR[3]=%e\n", RPAR[3]); |
---|
| 341 | return ros_ErrorMsg(-4,Tstart,ZERO); |
---|
| 342 | } /* end if */ |
---|
| 343 | if (RPAR[4] == ZERO) |
---|
| 344 | FacMax = (KPP_REAL)6.0; |
---|
| 345 | else |
---|
| 346 | FacMax = RPAR[4]; |
---|
| 347 | if (RPAR[4] < ZERO) { |
---|
| 348 | printf("\n User-selected FacMax: RPAR[4]=%e\n", RPAR[4]); |
---|
| 349 | return ros_ErrorMsg(-4,Tstart,ZERO); |
---|
| 350 | } /* end if */ |
---|
| 351 | /*~~~> FacRej: Factor to decrease step after 2 succesive rejections */ |
---|
| 352 | if (RPAR[5] == ZERO) |
---|
| 353 | FacRej = (KPP_REAL)0.1; |
---|
| 354 | else |
---|
| 355 | FacRej = RPAR[5]; |
---|
| 356 | if (RPAR[5] < ZERO) { |
---|
| 357 | printf("\n User-selected FacRej: RPAR[5]=%e\n", RPAR[5]); |
---|
| 358 | return ros_ErrorMsg(-4,Tstart,ZERO); |
---|
| 359 | } /* end if */ |
---|
| 360 | /*~~~> FacSafe: Safety Factor in the computation of new step size */ |
---|
| 361 | if (RPAR[6] == ZERO) |
---|
| 362 | FacSafe = (KPP_REAL)0.9; |
---|
| 363 | else |
---|
| 364 | FacSafe = RPAR[6]; |
---|
| 365 | if (RPAR[6] < ZERO) { |
---|
| 366 | printf("\n User-selected FacSafe: RPAR[6]=%e\n", RPAR[6]); |
---|
| 367 | return ros_ErrorMsg(-4,Tstart,ZERO); |
---|
| 368 | } /* end if */ |
---|
| 369 | /*~~~> Check if tolerances are reasonable */ |
---|
| 370 | for (i = 0; i < UplimTol; i++) { |
---|
| 371 | if ( (AbsTol[i] <= ZERO) || (RelTol[i] <= 10.0*Roundoff) |
---|
| 372 | || (RelTol[i] >= ONE) ) { |
---|
| 373 | printf("\n AbsTol[%d] = %e\n",i,AbsTol[i]); |
---|
| 374 | printf("\n RelTol[%d] = %e\n",i,RelTol[i]); |
---|
| 375 | return ros_ErrorMsg(-5,Tstart,ZERO); |
---|
| 376 | } /* end if */ |
---|
| 377 | } /* for */ |
---|
| 378 | |
---|
| 379 | |
---|
| 380 | /*~~~> Initialize the particular Rosenbrock method */ |
---|
| 381 | switch (Method) { |
---|
| 382 | case 1: |
---|
| 383 | Ros2(&ros_S, ros_A, ros_C, ros_M, ros_E, |
---|
| 384 | ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name); |
---|
| 385 | break; |
---|
| 386 | case 2: |
---|
| 387 | Ros3(&ros_S, ros_A, ros_C, ros_M, ros_E, |
---|
| 388 | ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name); |
---|
| 389 | break; |
---|
| 390 | case 3: |
---|
| 391 | Ros4(&ros_S, ros_A, ros_C, ros_M, ros_E, |
---|
| 392 | ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name); |
---|
| 393 | break; |
---|
| 394 | case 4: |
---|
| 395 | Rodas3(&ros_S, ros_A, ros_C, ros_M, ros_E, |
---|
| 396 | ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name); |
---|
| 397 | break; |
---|
| 398 | case 5: |
---|
| 399 | Rodas4(&ros_S, ros_A, ros_C, ros_M, ros_E, |
---|
| 400 | ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name); |
---|
| 401 | break; |
---|
| 402 | default: |
---|
| 403 | printf("\n Unknown Rosenbrock method: IPAR[3]= %d", Method); |
---|
| 404 | return ros_ErrorMsg(-2,Tstart,ZERO); |
---|
| 405 | } /* end switch */ |
---|
| 406 | |
---|
| 407 | /*~~~> Rosenbrock method */ |
---|
| 408 | IERR = RosenbrockIntegrator( Y,Tstart,Tend, |
---|
| 409 | AbsTol, RelTol, |
---|
| 410 | ode_Fun,ode_Jac , |
---|
| 411 | /* Rosenbrock method coefficients */ |
---|
| 412 | ros_S, ros_M, ros_E, ros_A, ros_C, |
---|
| 413 | ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, |
---|
| 414 | /* Integration parameters */ |
---|
| 415 | Autonomous, VectorTol, Max_no_steps, |
---|
| 416 | Roundoff, Hmin, Hmax, Hstart, |
---|
| 417 | FacMin, FacMax, FacRej, FacSafe, |
---|
| 418 | /* Output parameters */ |
---|
| 419 | &Texit, &Hexit ); |
---|
| 420 | |
---|
| 421 | |
---|
| 422 | /*~~~> Collect run statistics */ |
---|
| 423 | IPAR[10] = Nfun; |
---|
| 424 | IPAR[11] = Njac; |
---|
| 425 | IPAR[12] = Nstp; |
---|
| 426 | IPAR[13] = Nacc; |
---|
| 427 | IPAR[14] = Nrej; |
---|
| 428 | IPAR[15] = Ndec; |
---|
| 429 | IPAR[16] = Nsol; |
---|
| 430 | IPAR[17] = Nsng; |
---|
| 431 | /*~~~> Last T and H */ |
---|
| 432 | RPAR[10] = Texit; |
---|
| 433 | RPAR[11] = Hexit; |
---|
| 434 | |
---|
| 435 | return IERR; |
---|
| 436 | |
---|
| 437 | } /* Rosenbrock */ |
---|
| 438 | |
---|
| 439 | |
---|
| 440 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 441 | int RosenbrockIntegrator( |
---|
| 442 | /*~~~> Input: the initial condition at Tstart; Output: the solution at T */ |
---|
| 443 | KPP_REAL Y[], |
---|
| 444 | /*~~~> Input: integration interval */ |
---|
| 445 | KPP_REAL Tstart, KPP_REAL Tend , |
---|
| 446 | /*~~~> Input: tolerances */ |
---|
| 447 | KPP_REAL AbsTol[], KPP_REAL RelTol[], |
---|
| 448 | /*~~~> Input: ode function and its Jacobian */ |
---|
| 449 | void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), |
---|
| 450 | void (*ode_Jac)(KPP_REAL, KPP_REAL [], KPP_REAL []) , |
---|
| 451 | /*~~~> Input: The Rosenbrock method parameters */ |
---|
| 452 | int ros_S, |
---|
| 453 | KPP_REAL ros_M[], KPP_REAL ros_E[], |
---|
| 454 | KPP_REAL ros_A[], KPP_REAL ros_C[], |
---|
| 455 | KPP_REAL ros_Alpha[],KPP_REAL ros_Gamma[], |
---|
| 456 | KPP_REAL ros_ELO, char ros_NewF[], |
---|
| 457 | /*~~~> Input: integration parameters */ |
---|
| 458 | char Autonomous, char VectorTol, |
---|
| 459 | int Max_no_steps, |
---|
| 460 | KPP_REAL Roundoff, KPP_REAL Hmin, KPP_REAL Hmax, KPP_REAL Hstart, |
---|
| 461 | KPP_REAL FacMin, KPP_REAL FacMax, KPP_REAL FacRej, KPP_REAL FacSafe, |
---|
| 462 | /*~~~> Output: time at which the solution is returned (T=Tend if success) |
---|
| 463 | and last accepted step */ |
---|
| 464 | KPP_REAL *Texit, KPP_REAL *Hexit ) |
---|
| 465 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 466 | Template for the implementation of a generic Rosenbrock method |
---|
| 467 | defined by ros_S (no of stages) and coefficients ros_{A,C,M,E,Alpha,Gamma} |
---|
| 468 | |
---|
| 469 | returned value: IERR, indicator of success (if positive) |
---|
| 470 | or failure (if negative) |
---|
| 471 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 472 | { |
---|
| 473 | KPP_REAL Ynew[KPP_NVAR], Fcn0[KPP_NVAR], Fcn[KPP_NVAR], |
---|
| 474 | dFdT[KPP_NVAR], |
---|
| 475 | Jac0[KPP_LU_NONZERO], Ghimj[KPP_LU_NONZERO]; |
---|
| 476 | KPP_REAL K[KPP_NVAR*ros_S]; |
---|
| 477 | KPP_REAL H, T, Hnew, HC, HG, Fac, Tau; |
---|
| 478 | KPP_REAL Err, Yerr[KPP_NVAR]; |
---|
| 479 | int Pivot[KPP_NVAR], Direction, ioffset, j, istage; |
---|
| 480 | char RejectLastH, RejectMoreH; |
---|
| 481 | |
---|
| 482 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 483 | |
---|
| 484 | /*~~~> INITIAL PREPARATIONS */ |
---|
| 485 | T = Tstart; |
---|
| 486 | *Hexit = 0.0; |
---|
| 487 | H = MIN(Hstart,Hmax); |
---|
| 488 | if (ABS(H) <= 10.0*Roundoff) |
---|
| 489 | H = DeltaMin; |
---|
| 490 | |
---|
| 491 | if (Tend >= Tstart) { |
---|
| 492 | Direction = +1; |
---|
| 493 | } else { |
---|
| 494 | Direction = -1; |
---|
| 495 | } /* end if */ |
---|
| 496 | |
---|
| 497 | RejectLastH=0; RejectMoreH=0; |
---|
| 498 | |
---|
| 499 | /*~~~> Time loop begins below */ |
---|
| 500 | |
---|
| 501 | while ( ( (Direction > 0) && ((T-Tend)+Roundoff <= ZERO) ) |
---|
| 502 | || ( (Direction < 0) && ((Tend-T)+Roundoff <= ZERO) ) ) { |
---|
| 503 | |
---|
| 504 | if ( Nstp > Max_no_steps ) { /* Too many steps */ |
---|
| 505 | *Texit = T; |
---|
| 506 | return ros_ErrorMsg(-6,T,H); |
---|
| 507 | } |
---|
| 508 | if ( ((T+0.1*H) == T) || (H <= Roundoff) ) { /* Step size too small */ |
---|
| 509 | *Texit = T; |
---|
| 510 | return ros_ErrorMsg(-7,T,H); |
---|
| 511 | } |
---|
| 512 | |
---|
| 513 | /*~~~> Limit H if necessary to avoid going beyond Tend */ |
---|
| 514 | *Hexit = H; |
---|
| 515 | H = MIN(H,ABS(Tend-T)); |
---|
| 516 | |
---|
| 517 | /*~~~> Compute the function at current time */ |
---|
| 518 | (*ode_Fun)(T,Y,Fcn0); |
---|
| 519 | |
---|
| 520 | /*~~~> Compute the function derivative with respect to T */ |
---|
| 521 | if (!Autonomous) |
---|
| 522 | ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, ode_Fun, dFdT ); |
---|
| 523 | |
---|
| 524 | /*~~~> Compute the Jacobian at current time */ |
---|
| 525 | (*ode_Jac)(T,Y,Jac0); |
---|
| 526 | |
---|
| 527 | /*~~~> Repeat step calculation until current step accepted */ |
---|
| 528 | while (1) { /* WHILE STEP NOT ACCEPTED */ |
---|
| 529 | |
---|
| 530 | |
---|
| 531 | if( ros_PrepareMatrix( &H, Direction, ros_Gamma[0], |
---|
| 532 | Jac0, Ghimj, Pivot) ) { /* More than 5 consecutive failed decompositions */ |
---|
| 533 | *Texit = T; |
---|
| 534 | return ros_ErrorMsg(-8,T,H); |
---|
| 535 | } |
---|
| 536 | |
---|
| 537 | /*~~~> Compute the stages */ |
---|
| 538 | for (istage = 1; istage <= ros_S; istage++) { |
---|
| 539 | |
---|
| 540 | /* Current istage offset. Current istage vector is K[ioffset:ioffset+KPP_NVAR-1] */ |
---|
| 541 | ioffset = KPP_NVAR*(istage-1); |
---|
| 542 | |
---|
| 543 | /* For the 1st istage the function has been computed previously */ |
---|
| 544 | if ( istage == 1 ) |
---|
| 545 | WCOPY(KPP_NVAR,Fcn0,1,Fcn,1); |
---|
| 546 | else { /* istage>1 and a new function evaluation is needed at current istage */ |
---|
| 547 | if ( ros_NewF[istage-1] ) { |
---|
| 548 | WCOPY(KPP_NVAR,Y,1,Ynew,1); |
---|
| 549 | for (j = 1; j <= istage-1; j++) |
---|
| 550 | WAXPY(KPP_NVAR,ros_A[(istage-1)*(istage-2)/2+j-1], |
---|
| 551 | &K[KPP_NVAR*(j-1)],1,Ynew,1); |
---|
| 552 | Tau = T + ros_Alpha[istage-1]*Direction*H; |
---|
| 553 | (*ode_Fun)(Tau,Ynew,Fcn); |
---|
| 554 | } /*end if ros_NewF(istage)*/ |
---|
| 555 | } /* end if istage */ |
---|
| 556 | |
---|
| 557 | WCOPY(KPP_NVAR,Fcn,1,&K[ioffset],1); |
---|
| 558 | for (j = 1; j <= istage-1; j++) { |
---|
| 559 | HC = ros_C[(istage-1)*(istage-2)/2+j-1]/(Direction*H); |
---|
| 560 | WAXPY(KPP_NVAR,HC,&K[KPP_NVAR*(j-1)],1,&K[ioffset],1); |
---|
| 561 | } /* for j */ |
---|
| 562 | |
---|
| 563 | if ((!Autonomous) && (ros_Gamma[istage-1])) { |
---|
| 564 | HG = Direction*H*ros_Gamma[istage-1]; |
---|
| 565 | WAXPY(KPP_NVAR,HG,dFdT,1,&K[ioffset],1); |
---|
| 566 | } /* end if !Autonomous */ |
---|
| 567 | |
---|
| 568 | SolveTemplate(Ghimj, Pivot, &K[ioffset]); |
---|
| 569 | |
---|
| 570 | } /* for istage */ |
---|
| 571 | |
---|
| 572 | |
---|
| 573 | /*~~~> Compute the new solution */ |
---|
| 574 | WCOPY(KPP_NVAR,Y,1,Ynew,1); |
---|
| 575 | for (j=1; j<=ros_S; j++) |
---|
| 576 | WAXPY(KPP_NVAR,ros_M[j-1],&K[KPP_NVAR*(j-1)],1,Ynew,1); |
---|
| 577 | |
---|
| 578 | /*~~~> Compute the error estimation */ |
---|
| 579 | WSCAL(KPP_NVAR,ZERO,Yerr,1); |
---|
| 580 | for (j=1; j<=ros_S; j++) |
---|
| 581 | WAXPY(KPP_NVAR,ros_E[j-1],&K[KPP_NVAR*(j-1)],1,Yerr,1); |
---|
| 582 | Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ); |
---|
| 583 | |
---|
| 584 | /*~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax */ |
---|
| 585 | Fac = MIN(FacMax,MAX(FacMin,FacSafe/pow(Err,ONE/ros_ELO))); |
---|
| 586 | Hnew = H*Fac; |
---|
| 587 | |
---|
| 588 | /*~~~> Check the error magnitude and adjust step size */ |
---|
| 589 | Nstp++; |
---|
| 590 | if ( (Err <= ONE) || (H <= Hmin) ) { /*~~~> Accept step */ |
---|
| 591 | Nacc++; |
---|
| 592 | WCOPY(KPP_NVAR,Ynew,1,Y,1); |
---|
| 593 | T += Direction*H; |
---|
| 594 | Hnew = MAX(Hmin,MIN(Hnew,Hmax)); |
---|
| 595 | /* No step size increase after a rejected step */ |
---|
| 596 | if (RejectLastH) |
---|
| 597 | Hnew = MIN(Hnew,H); |
---|
| 598 | RejectLastH = 0; RejectMoreH = 0; |
---|
| 599 | H = Hnew; |
---|
| 600 | break; /* EXIT THE LOOP: WHILE STEP NOT ACCEPTED */ |
---|
| 601 | } else { /*~~~> Reject step */ |
---|
| 602 | if (Nacc >= 1) |
---|
| 603 | Nrej++; |
---|
| 604 | if (RejectMoreH) |
---|
| 605 | Hnew=H*FacRej; |
---|
| 606 | RejectMoreH = RejectLastH; RejectLastH = 1; |
---|
| 607 | H = Hnew; |
---|
| 608 | } /* end if Err <= 1 */ |
---|
| 609 | |
---|
| 610 | } /* while LOOP: WHILE STEP NOT ACCEPTED */ |
---|
| 611 | |
---|
| 612 | } /* while: time loop */ |
---|
| 613 | |
---|
| 614 | /*~~~> The integration was successful */ |
---|
| 615 | *Texit = T; |
---|
| 616 | return 1; |
---|
| 617 | |
---|
| 618 | } /* RosenbrockIntegrator */ |
---|
| 619 | |
---|
| 620 | |
---|
| 621 | |
---|
| 622 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 623 | KPP_REAL ros_ErrorNorm ( |
---|
| 624 | /*~~~> Input arguments */ |
---|
| 625 | KPP_REAL Y[], KPP_REAL Ynew[], KPP_REAL Yerr[], |
---|
| 626 | KPP_REAL AbsTol[], KPP_REAL RelTol[], |
---|
| 627 | char VectorTol ) |
---|
| 628 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 629 | Computes and returns the "scaled norm" of the error vector Yerr |
---|
| 630 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 631 | { |
---|
| 632 | /*~~~> Local variables */ |
---|
| 633 | KPP_REAL Err, Scale, Ymax; |
---|
| 634 | int i; |
---|
| 635 | |
---|
| 636 | Err = ZERO; |
---|
| 637 | for (i=0; i<KPP_NVAR; i++) { |
---|
| 638 | Ymax = MAX(ABS(Y[i]),ABS(Ynew[i])); |
---|
| 639 | if (VectorTol) { |
---|
| 640 | Scale = AbsTol[i]+RelTol[i]*Ymax; |
---|
| 641 | } else { |
---|
| 642 | Scale = AbsTol[0]+RelTol[0]*Ymax; |
---|
| 643 | } /* end if */ |
---|
| 644 | Err = Err+(Yerr[i]*Yerr[i])/(Scale*Scale); |
---|
| 645 | } /* for i */ |
---|
| 646 | Err = SQRT(Err/(KPP_REAL)KPP_NVAR); |
---|
| 647 | |
---|
| 648 | return Err; |
---|
| 649 | |
---|
| 650 | } /* ros_ErrorNorm */ |
---|
| 651 | |
---|
| 652 | |
---|
| 653 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 654 | void ros_FunTimeDerivative ( |
---|
| 655 | /*~~~> Input arguments: */ |
---|
| 656 | KPP_REAL T, KPP_REAL Roundoff, |
---|
| 657 | KPP_REAL Y[], KPP_REAL Fcn0[], |
---|
| 658 | void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), |
---|
| 659 | /*~~~> Output arguments: */ |
---|
| 660 | KPP_REAL dFdT[] ) |
---|
| 661 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 662 | The time partial derivative of the function by finite differences |
---|
| 663 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 664 | { |
---|
| 665 | /*~~~> Local variables */ |
---|
| 666 | KPP_REAL Delta; |
---|
| 667 | |
---|
| 668 | Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)); |
---|
| 669 | (*ode_Fun)(T+Delta,Y,dFdT); |
---|
| 670 | WAXPY(KPP_NVAR,(-ONE),Fcn0,1,dFdT,1); |
---|
| 671 | WSCAL(KPP_NVAR,(ONE/Delta),dFdT,1); |
---|
| 672 | |
---|
| 673 | } /* ros_FunTimeDerivative */ |
---|
| 674 | |
---|
| 675 | |
---|
| 676 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 677 | char ros_PrepareMatrix ( |
---|
| 678 | /* Inout argument: (step size is decreased when LU fails) */ |
---|
| 679 | KPP_REAL* H, |
---|
| 680 | /* Input arguments: */ |
---|
| 681 | int Direction, KPP_REAL gam, KPP_REAL Jac0[], |
---|
| 682 | /* Output arguments: */ |
---|
| 683 | KPP_REAL Ghimj[], int Pivot[] ) |
---|
| 684 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 685 | Prepares the LHS matrix for stage calculations |
---|
| 686 | 1. Construct Ghimj = 1/(H*ham) - Jac0 |
---|
| 687 | "(Gamma H) Inverse Minus Jacobian" |
---|
| 688 | 2. Repeat LU decomposition of Ghimj until successful. |
---|
| 689 | -half the step size if LU decomposition fails and retry |
---|
| 690 | -exit after 5 consecutive fails |
---|
| 691 | |
---|
| 692 | Return value: Singular (true=1=failed_LU or false=0=successful_LU) |
---|
| 693 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 694 | { |
---|
| 695 | /*~~~> Local variables */ |
---|
| 696 | int i, ising, Nconsecutive; |
---|
| 697 | KPP_REAL ghinv; |
---|
| 698 | |
---|
| 699 | Nconsecutive = 0; |
---|
| 700 | |
---|
| 701 | while (1) { /* while Singular */ |
---|
| 702 | |
---|
| 703 | /*~~~> Construct Ghimj = 1/(H*ham) - Jac0 */ |
---|
| 704 | WCOPY(KPP_LU_NONZERO,Jac0,1,Ghimj,1); |
---|
| 705 | WSCAL(KPP_LU_NONZERO,(-ONE),Ghimj,1); |
---|
| 706 | ghinv = ONE/(Direction*(*H)*gam); |
---|
| 707 | for (i=0; i<KPP_NVAR; i++) { |
---|
| 708 | Ghimj[LU_DIAG[i]] = Ghimj[LU_DIAG[i]]+ghinv; |
---|
| 709 | } /* for i */ |
---|
| 710 | /*~~~> Compute LU decomposition */ |
---|
| 711 | DecompTemplate( Ghimj, Pivot, &ising ); |
---|
| 712 | if (ising == 0) { |
---|
| 713 | /*~~~> if successful done */ |
---|
| 714 | return 0; /* Singular = false */ |
---|
| 715 | } else { /* ising .ne. 0 */ |
---|
| 716 | /*~~~> if unsuccessful half the step size; if 5 consecutive fails return */ |
---|
| 717 | Nsng++; Nconsecutive++; |
---|
| 718 | printf("\nWarning: LU Decomposition returned ising = %d\n",ising); |
---|
| 719 | if (Nconsecutive <= 5) { /* Less than 5 consecutive failed LUs */ |
---|
| 720 | *H = (*H)*HALF; |
---|
| 721 | } else { /* More than 5 consecutive failed LUs */ |
---|
| 722 | return 1; /* Singular = true */ |
---|
| 723 | } /* end if Nconsecutive */ |
---|
| 724 | } /* end if ising */ |
---|
| 725 | |
---|
| 726 | } /* while Singular */ |
---|
| 727 | |
---|
| 728 | } /* ros_PrepareMatrix */ |
---|
| 729 | |
---|
| 730 | |
---|
| 731 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 732 | int ros_ErrorMsg(int Code, KPP_REAL T, KPP_REAL H) |
---|
| 733 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 734 | Handles all error messages and returns IERR = error Code |
---|
| 735 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 736 | { |
---|
| 737 | printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"); |
---|
| 738 | printf("\nForced exit from Rosenbrock due to the following error:\n"); |
---|
| 739 | |
---|
| 740 | switch (Code) { |
---|
| 741 | case -1: |
---|
| 742 | printf("--> Improper value for maximal no of steps"); break; |
---|
| 743 | case -2: |
---|
| 744 | printf("--> Selected Rosenbrock method not implemented"); break; |
---|
| 745 | case -3: |
---|
| 746 | printf("--> Hmin/Hmax/Hstart must be positive"); break; |
---|
| 747 | case -4: |
---|
| 748 | printf("--> FacMin/FacMax/FacRej must be positive"); break; |
---|
| 749 | case -5: |
---|
| 750 | printf("--> Improper tolerance values"); break; |
---|
| 751 | case -6: |
---|
| 752 | printf("--> No of steps exceeds maximum bound"); break; |
---|
| 753 | case -7: |
---|
| 754 | printf("--> Step size too small (T + H/10 = T) or H < Roundoff"); break; |
---|
| 755 | case -8: |
---|
| 756 | printf("--> Matrix is repeatedly singular"); break; |
---|
| 757 | default: |
---|
| 758 | printf("Unknown Error code: %d ",Code); |
---|
| 759 | } /* end switch */ |
---|
| 760 | |
---|
| 761 | printf("\n Time = %15.7e, H = %15.7e",T,H); |
---|
| 762 | printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"); |
---|
| 763 | |
---|
| 764 | return Code; |
---|
| 765 | |
---|
| 766 | } /* ros_ErrorMsg */ |
---|
| 767 | |
---|
| 768 | |
---|
| 769 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 770 | void Ros2 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], |
---|
| 771 | KPP_REAL ros_M[], KPP_REAL ros_E[], |
---|
| 772 | KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], |
---|
| 773 | char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ) |
---|
| 774 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 775 | AN L-STABLE METHOD, 2 stages, order 2 |
---|
| 776 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 777 | { |
---|
| 778 | double g = (KPP_REAL)1.70710678118655; /* 1.0 + 1.0/SQRT(2.0) */ |
---|
| 779 | |
---|
| 780 | /*~~~> Name of the method */ |
---|
| 781 | strcpy(ros_Name, "ROS-2"); |
---|
| 782 | |
---|
| 783 | /*~~~> Number of stages */ |
---|
| 784 | *ros_S = 2; |
---|
| 785 | |
---|
| 786 | /*~~~> The coefficient matrices A and C are strictly lower triangular. |
---|
| 787 | The lower triangular (subdiagonal) elements are stored in row-wise order: |
---|
| 788 | A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc. |
---|
| 789 | The general mapping formula is: |
---|
| 790 | A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ] */ |
---|
| 791 | ros_A[0] = 1.0/g; |
---|
| 792 | |
---|
| 793 | /*~~~> C_{i,j} = ros_C[ (i-1)*(i-2)/2 + j -1] */ |
---|
| 794 | ros_C[0] = (-2.0)/g; |
---|
| 795 | |
---|
| 796 | /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE) |
---|
| 797 | or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */ |
---|
| 798 | ros_NewF[0] = 1; |
---|
| 799 | ros_NewF[1] = 1; |
---|
| 800 | |
---|
| 801 | /*~~~> M_i = Coefficients for new step solution */ |
---|
| 802 | ros_M[0]= (3.0)/(2.0*g); |
---|
| 803 | ros_M[1]= (1.0)/(2.0*g); |
---|
| 804 | |
---|
| 805 | /*~~~> E_i = Coefficients for error estimator */ |
---|
| 806 | ros_E[0] = 1.0/(2.0*g); |
---|
| 807 | ros_E[1] = 1.0/(2.0*g); |
---|
| 808 | |
---|
| 809 | /*~~~> ros_ELO = estimator of local order - the minimum between the |
---|
| 810 | ! main and the embedded scheme orders plus one */ |
---|
| 811 | *ros_ELO = (KPP_REAL)2.0; |
---|
| 812 | |
---|
| 813 | /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ |
---|
| 814 | ros_Alpha[0] = (KPP_REAL)0.0; |
---|
| 815 | ros_Alpha[1] = (KPP_REAL)1.0; |
---|
| 816 | |
---|
| 817 | /*~~~> Gamma_i = \sum_j gamma_{i,j} */ |
---|
| 818 | ros_Gamma[0] = g; |
---|
| 819 | ros_Gamma[1] = -g; |
---|
| 820 | |
---|
| 821 | } /* Ros2 */ |
---|
| 822 | |
---|
| 823 | |
---|
| 824 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 825 | void Ros3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], |
---|
| 826 | KPP_REAL ros_M[], KPP_REAL ros_E[], |
---|
| 827 | KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], |
---|
| 828 | char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ) |
---|
| 829 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 830 | AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations |
---|
| 831 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 832 | { |
---|
| 833 | /*~~~> Name of the method */ |
---|
| 834 | strcpy(ros_Name, "ROS-3"); |
---|
| 835 | |
---|
| 836 | /*~~~> Number of stages */ |
---|
| 837 | *ros_S = 3; |
---|
| 838 | |
---|
| 839 | /*~~~> The coefficient matrices A and C are strictly lower triangular. |
---|
| 840 | The lower triangular (subdiagonal) elements are stored in row-wise order: |
---|
| 841 | A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc. |
---|
| 842 | The general mapping formula is: |
---|
| 843 | A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ] */ |
---|
| 844 | ros_A[0]= (KPP_REAL)1.0; |
---|
| 845 | ros_A[1]= (KPP_REAL)1.0; |
---|
| 846 | ros_A[2]= (KPP_REAL)0.0; |
---|
| 847 | |
---|
| 848 | /*~~~> C_{i,j} = ros_C[ (i-1)*(i-2)/2 + j -1] */ |
---|
| 849 | ros_C[0] = (KPP_REAL)(-1.0156171083877702091975600115545); |
---|
| 850 | ros_C[1] = (KPP_REAL)4.0759956452537699824805835358067; |
---|
| 851 | ros_C[2] = (KPP_REAL)9.2076794298330791242156818474003; |
---|
| 852 | |
---|
| 853 | /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE) |
---|
| 854 | or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */ |
---|
| 855 | ros_NewF[0] = 1; |
---|
| 856 | ros_NewF[1] = 1; |
---|
| 857 | ros_NewF[2] = 0; |
---|
| 858 | |
---|
| 859 | /*~~~> M_i = Coefficients for new step solution */ |
---|
| 860 | ros_M[0] = (KPP_REAL)1.0; |
---|
| 861 | ros_M[1] = (KPP_REAL)6.1697947043828245592553615689730; |
---|
| 862 | ros_M[2] = (KPP_REAL)(-0.4277225654321857332623837380651); |
---|
| 863 | |
---|
| 864 | /*~~~> E_i = Coefficients for error estimator */ |
---|
| 865 | ros_E[0] = (KPP_REAL)0.5; |
---|
| 866 | ros_E[1] = (KPP_REAL)(-2.9079558716805469821718236208017); |
---|
| 867 | ros_E[2] = (KPP_REAL)0.2235406989781156962736090927619; |
---|
| 868 | |
---|
| 869 | /*~~~> ros_ELO = estimator of local order - the minimum between the |
---|
| 870 | ! main and the embedded scheme orders plus 1 */ |
---|
| 871 | *ros_ELO = (KPP_REAL)3.0; |
---|
| 872 | |
---|
| 873 | /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ |
---|
| 874 | ros_Alpha[0]= (KPP_REAL)0.0; |
---|
| 875 | ros_Alpha[1]= (KPP_REAL)0.43586652150845899941601945119356; |
---|
| 876 | ros_Alpha[2]= (KPP_REAL)0.43586652150845899941601945119356; |
---|
| 877 | |
---|
| 878 | /*~~~> Gamma_i = \sum_j gamma_{i,j} */ |
---|
| 879 | ros_Gamma[0]= (KPP_REAL)0.43586652150845899941601945119356; |
---|
| 880 | ros_Gamma[1]= (KPP_REAL)0.24291996454816804366592249683314; |
---|
| 881 | ros_Gamma[2]= (KPP_REAL)2.1851380027664058511513169485832; |
---|
| 882 | |
---|
| 883 | } /* Ros3 */ |
---|
| 884 | |
---|
| 885 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 886 | |
---|
| 887 | |
---|
| 888 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 889 | void Ros4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], |
---|
| 890 | KPP_REAL ros_M[], KPP_REAL ros_E[], |
---|
| 891 | KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], |
---|
| 892 | char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ) |
---|
| 893 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 894 | L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES |
---|
| 895 | L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3 |
---|
| 896 | |
---|
| 897 | E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL |
---|
| 898 | EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. |
---|
| 899 | SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, |
---|
| 900 | SPRINGER-VERLAG (1990) |
---|
| 901 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 902 | { |
---|
| 903 | /*~~~> Name of the method */ |
---|
| 904 | strcpy(ros_Name, "ROS-4"); |
---|
| 905 | |
---|
| 906 | /*~~~> Number of stages */ |
---|
| 907 | *ros_S = 4; |
---|
| 908 | |
---|
| 909 | /*~~~> The coefficient matrices A and C are strictly lower triangular. |
---|
| 910 | The lower triangular (subdiagonal) elements are stored in row-wise order: |
---|
| 911 | A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc. |
---|
| 912 | The general mapping formula is: |
---|
| 913 | A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ] */ |
---|
| 914 | ros_A[0] = (KPP_REAL)0.2000000000000000e+01; |
---|
| 915 | ros_A[1] = (KPP_REAL)0.1867943637803922e+01; |
---|
| 916 | ros_A[2] = (KPP_REAL)0.2344449711399156; |
---|
| 917 | ros_A[3] = ros_A[1]; |
---|
| 918 | ros_A[4] = ros_A[2]; |
---|
| 919 | ros_A[5] = (KPP_REAL)0.0; |
---|
| 920 | |
---|
| 921 | /*~~~> C(i,j) = (KPP_REAL)ros_C( (i-1)*(i-2)/2 + j ) */ |
---|
| 922 | ros_C[0] = (KPP_REAL)(-0.7137615036412310e+01); |
---|
| 923 | ros_C[1] = (KPP_REAL)( 0.2580708087951457e+01); |
---|
| 924 | ros_C[2] = (KPP_REAL)( 0.6515950076447975); |
---|
| 925 | ros_C[3] = (KPP_REAL)(-0.2137148994382534e+01); |
---|
| 926 | ros_C[4] = (KPP_REAL)(-0.3214669691237626); |
---|
| 927 | ros_C[5] = (KPP_REAL)(-0.6949742501781779); |
---|
| 928 | |
---|
| 929 | /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE) |
---|
| 930 | or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */ |
---|
| 931 | ros_NewF[0] = 1; |
---|
| 932 | ros_NewF[1] = 1; |
---|
| 933 | ros_NewF[2] = 1; |
---|
| 934 | ros_NewF[3] = 0; |
---|
| 935 | |
---|
| 936 | /*~~~> M_i = Coefficients for new step solution */ |
---|
| 937 | ros_M[0] = (KPP_REAL)0.2255570073418735e+01; |
---|
| 938 | ros_M[1] = (KPP_REAL)0.2870493262186792; |
---|
| 939 | ros_M[2] = (KPP_REAL)0.4353179431840180; |
---|
| 940 | ros_M[3] = (KPP_REAL)0.1093502252409163e+01; |
---|
| 941 | |
---|
| 942 | /*~~~> E_i = Coefficients for error estimator */ |
---|
| 943 | ros_E[0] = (KPP_REAL)(-0.2815431932141155); |
---|
| 944 | ros_E[1] = (KPP_REAL)(-0.7276199124938920e-01); |
---|
| 945 | ros_E[2] = (KPP_REAL)(-0.1082196201495311); |
---|
| 946 | ros_E[3] = (KPP_REAL)(-0.1093502252409163e+01); |
---|
| 947 | |
---|
| 948 | /*~~~> ros_ELO = estimator of local order - the minimum between the |
---|
| 949 | ! main and the embedded scheme orders plus 1 */ |
---|
| 950 | *ros_ELO = (KPP_REAL)4.0; |
---|
| 951 | |
---|
| 952 | /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ |
---|
| 953 | ros_Alpha[0] = (KPP_REAL)0.0; |
---|
| 954 | ros_Alpha[1] = (KPP_REAL)0.1145640000000000e+01; |
---|
| 955 | ros_Alpha[2] = (KPP_REAL)0.6552168638155900; |
---|
| 956 | ros_Alpha[3] = (KPP_REAL)ros_Alpha[2]; |
---|
| 957 | |
---|
| 958 | /*~~~> Gamma_i = \sum_j gamma_{i,j} */ |
---|
| 959 | ros_Gamma[0] = (KPP_REAL)( 0.5728200000000000); |
---|
| 960 | ros_Gamma[1] = (KPP_REAL)(-0.1769193891319233e+01); |
---|
| 961 | ros_Gamma[2] = (KPP_REAL)( 0.7592633437920482); |
---|
| 962 | ros_Gamma[3] = (KPP_REAL)(-0.1049021087100450); |
---|
| 963 | |
---|
| 964 | } /* Ros4 */ |
---|
| 965 | |
---|
| 966 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 967 | void Rodas3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], |
---|
| 968 | KPP_REAL ros_M[], KPP_REAL ros_E[], |
---|
| 969 | KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], |
---|
| 970 | char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ) |
---|
| 971 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 972 | --- A STIFFLY-STABLE METHOD, 4 stages, order 3 |
---|
| 973 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 974 | { |
---|
| 975 | /*~~~> Name of the method */ |
---|
| 976 | strcpy(ros_Name, "RODAS-3"); |
---|
| 977 | |
---|
| 978 | /*~~~> Number of stages */ |
---|
| 979 | *ros_S = 4; |
---|
| 980 | |
---|
| 981 | /*~~~> The coefficient matrices A and C are strictly lower triangular. |
---|
| 982 | The lower triangular (subdiagonal) elements are stored in row-wise order: |
---|
| 983 | A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc. |
---|
| 984 | The general mapping formula is: |
---|
| 985 | A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ] */ |
---|
| 986 | ros_A[0] = (KPP_REAL)0.0; |
---|
| 987 | ros_A[1] = (KPP_REAL)2.0; |
---|
| 988 | ros_A[2] = (KPP_REAL)0.0; |
---|
| 989 | ros_A[3] = (KPP_REAL)2.0; |
---|
| 990 | ros_A[4] = (KPP_REAL)0.0; |
---|
| 991 | ros_A[5] = (KPP_REAL)1.0; |
---|
| 992 | |
---|
| 993 | /*~~~> C_{i,j} = ros_C[ (i-1)*(i-2)/2 + j -1] */ |
---|
| 994 | ros_C[0] = (KPP_REAL)4.0; |
---|
| 995 | ros_C[1] = (KPP_REAL)1.0; |
---|
| 996 | ros_C[2] = (KPP_REAL)(-1.0); |
---|
| 997 | ros_C[3] = (KPP_REAL)1.0; |
---|
| 998 | ros_C[4] = (KPP_REAL)(-1.0); |
---|
| 999 | ros_C[5] = (KPP_REAL)(-2.66666666666667); /* -8/3 */ |
---|
| 1000 | |
---|
| 1001 | /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE) |
---|
| 1002 | or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */ |
---|
| 1003 | ros_NewF[0] = 1; |
---|
| 1004 | ros_NewF[1] = 0; |
---|
| 1005 | ros_NewF[2] = 1; |
---|
| 1006 | ros_NewF[3] = 1; |
---|
| 1007 | |
---|
| 1008 | /*~~~> M_i = Coefficients for new step solution */ |
---|
| 1009 | ros_M[0] = (KPP_REAL)2.0; |
---|
| 1010 | ros_M[1] = (KPP_REAL)0.0; |
---|
| 1011 | ros_M[2] = (KPP_REAL)1.0; |
---|
| 1012 | ros_M[3] = (KPP_REAL)1.0; |
---|
| 1013 | |
---|
| 1014 | /*~~~> E_i = Coefficients for error estimator */ |
---|
| 1015 | ros_E[0] = (KPP_REAL)0.0; |
---|
| 1016 | ros_E[1] = (KPP_REAL)0.0; |
---|
| 1017 | ros_E[2] = (KPP_REAL)0.0; |
---|
| 1018 | ros_E[3] = (KPP_REAL)1.0; |
---|
| 1019 | |
---|
| 1020 | /*~~~> ros_ELO = estimator of local order - the minimum between the |
---|
| 1021 | ! main and the embedded scheme orders plus 1 */ |
---|
| 1022 | *ros_ELO = (KPP_REAL)3.0; |
---|
| 1023 | |
---|
| 1024 | /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ |
---|
| 1025 | ros_Alpha[0] = (KPP_REAL)0.0; |
---|
| 1026 | ros_Alpha[1] = (KPP_REAL)0.0; |
---|
| 1027 | ros_Alpha[2] = (KPP_REAL)1.0; |
---|
| 1028 | ros_Alpha[3] = (KPP_REAL)1.0; |
---|
| 1029 | |
---|
| 1030 | /*~~~> Gamma_i = \sum_j gamma_{i,j} */ |
---|
| 1031 | ros_Gamma[0] = (KPP_REAL)0.5; |
---|
| 1032 | ros_Gamma[1] = (KPP_REAL)1.5; |
---|
| 1033 | ros_Gamma[2] = (KPP_REAL)0.0; |
---|
| 1034 | ros_Gamma[3] = (KPP_REAL)0.0; |
---|
| 1035 | |
---|
| 1036 | } /* Rodas3 */ |
---|
| 1037 | |
---|
| 1038 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1039 | void Rodas4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], |
---|
| 1040 | KPP_REAL ros_M[], KPP_REAL ros_E[], |
---|
| 1041 | KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], |
---|
| 1042 | char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ) |
---|
| 1043 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1044 | STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES |
---|
| 1045 | |
---|
| 1046 | E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL |
---|
| 1047 | EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. |
---|
| 1048 | SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, |
---|
| 1049 | SPRINGER-VERLAG (1996) |
---|
| 1050 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1051 | { |
---|
| 1052 | /*~~~> Name of the method */ |
---|
| 1053 | strcpy(ros_Name, "RODAS-4"); |
---|
| 1054 | |
---|
| 1055 | /*~~~> Number of stages */ |
---|
| 1056 | *ros_S = 6; |
---|
| 1057 | |
---|
| 1058 | /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ |
---|
| 1059 | ros_Alpha[0] = (KPP_REAL)0.000; |
---|
| 1060 | ros_Alpha[1] = (KPP_REAL)0.386; |
---|
| 1061 | ros_Alpha[2] = (KPP_REAL)0.210; |
---|
| 1062 | ros_Alpha[3] = (KPP_REAL)0.630; |
---|
| 1063 | ros_Alpha[4] = (KPP_REAL)1.000; |
---|
| 1064 | ros_Alpha[5] = (KPP_REAL)1.000; |
---|
| 1065 | |
---|
| 1066 | /*~~~> Gamma_i = \sum_j gamma_{i,j} */ |
---|
| 1067 | ros_Gamma[0] = (KPP_REAL)0.2500000000000000; |
---|
| 1068 | ros_Gamma[1] = (KPP_REAL)(-0.1043000000000000); |
---|
| 1069 | ros_Gamma[2] = (KPP_REAL)0.1035000000000000; |
---|
| 1070 | ros_Gamma[3] = (KPP_REAL)(-0.3620000000000023e-01); |
---|
| 1071 | ros_Gamma[4] = (KPP_REAL)0.0; |
---|
| 1072 | ros_Gamma[5] = (KPP_REAL)0.0; |
---|
| 1073 | |
---|
| 1074 | /*~~~> The coefficient matrices A and C are strictly lower triangular. |
---|
| 1075 | The lower triangular (subdiagonal) elements are stored in row-wise order: |
---|
| 1076 | A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc. |
---|
| 1077 | The general mapping formula is: A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ] */ |
---|
| 1078 | ros_A[0] = (KPP_REAL)0.1544000000000000e+01; |
---|
| 1079 | ros_A[1] = (KPP_REAL)0.9466785280815826; |
---|
| 1080 | ros_A[2] = (KPP_REAL)0.2557011698983284; |
---|
| 1081 | ros_A[3] = (KPP_REAL)0.3314825187068521e+01; |
---|
| 1082 | ros_A[4] = (KPP_REAL)0.2896124015972201e+01; |
---|
| 1083 | ros_A[5] = (KPP_REAL)0.9986419139977817; |
---|
| 1084 | ros_A[6] = (KPP_REAL)0.1221224509226641e+01; |
---|
| 1085 | ros_A[7] = (KPP_REAL)0.6019134481288629e+01; |
---|
| 1086 | ros_A[8] = (KPP_REAL)0.1253708332932087e+02; |
---|
| 1087 | ros_A[9] = (KPP_REAL)(-0.6878860361058950); |
---|
| 1088 | ros_A[10] = ros_A[6]; |
---|
| 1089 | ros_A[11] = ros_A[7]; |
---|
| 1090 | ros_A[12] = ros_A[8]; |
---|
| 1091 | ros_A[13] = ros_A[9]; |
---|
| 1092 | ros_A[14] = (KPP_REAL)1.0; |
---|
| 1093 | |
---|
| 1094 | /*~~~> C_{i,j} = ros_C[ (i-1)*(i-2)/2 + j -1] */ |
---|
| 1095 | ros_C[0] = (KPP_REAL)(-0.5668800000000000e+01); |
---|
| 1096 | ros_C[1] = (KPP_REAL)(-0.2430093356833875e+01); |
---|
| 1097 | ros_C[2] = (KPP_REAL)(-0.2063599157091915); |
---|
| 1098 | ros_C[3] = (KPP_REAL)(-0.1073529058151375); |
---|
| 1099 | ros_C[4] = (KPP_REAL)(-0.9594562251023355e+01); |
---|
| 1100 | ros_C[5] = (KPP_REAL)(-0.2047028614809616e+02); |
---|
| 1101 | ros_C[6] = (KPP_REAL)( 0.7496443313967647e+01); |
---|
| 1102 | ros_C[7] = (KPP_REAL)(-0.1024680431464352e+02); |
---|
| 1103 | ros_C[8] = (KPP_REAL)(-0.3399990352819905e+02); |
---|
| 1104 | ros_C[9] = (KPP_REAL)( 0.1170890893206160e+02); |
---|
| 1105 | ros_C[10] = (KPP_REAL)( 0.8083246795921522e+01); |
---|
| 1106 | ros_C[11] = (KPP_REAL)(-0.7981132988064893e+01); |
---|
| 1107 | ros_C[12] = (KPP_REAL)(-0.3152159432874371e+02); |
---|
| 1108 | ros_C[13] = (KPP_REAL)( 0.1631930543123136e+02); |
---|
| 1109 | ros_C[14] = (KPP_REAL)(-0.6058818238834054e+01); |
---|
| 1110 | |
---|
| 1111 | /*~~~> M_i = Coefficients for new step solution */ |
---|
| 1112 | ros_M[0] = ros_A[6]; |
---|
| 1113 | ros_M[1] = ros_A[7]; |
---|
| 1114 | ros_M[2] = ros_A[8]; |
---|
| 1115 | ros_M[3] = ros_A[9]; |
---|
| 1116 | ros_M[4] = (KPP_REAL)1.0; |
---|
| 1117 | ros_M[5] = (KPP_REAL)1.0; |
---|
| 1118 | |
---|
| 1119 | /*~~~> E_i = Coefficients for error estimator */ |
---|
| 1120 | ros_E[0] = (KPP_REAL)0.0; |
---|
| 1121 | ros_E[1] = (KPP_REAL)0.0; |
---|
| 1122 | ros_E[2] = (KPP_REAL)0.0; |
---|
| 1123 | ros_E[3] = (KPP_REAL)0.0; |
---|
| 1124 | ros_E[4] = (KPP_REAL)0.0; |
---|
| 1125 | ros_E[5] = (KPP_REAL)1.0; |
---|
| 1126 | |
---|
| 1127 | /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE) |
---|
| 1128 | or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */ |
---|
| 1129 | ros_NewF[0] = 1; |
---|
| 1130 | ros_NewF[1] = 1; |
---|
| 1131 | ros_NewF[2] = 1; |
---|
| 1132 | ros_NewF[3] = 1; |
---|
| 1133 | ros_NewF[4] = 1; |
---|
| 1134 | ros_NewF[5] = 1; |
---|
| 1135 | |
---|
| 1136 | /*~~~> ros_ELO = estimator of local order - the minimum between the |
---|
| 1137 | ! main and the embedded scheme orders plus 1 */ |
---|
| 1138 | *ros_ELO = (KPP_REAL)4.0; |
---|
| 1139 | |
---|
| 1140 | } /* Rodas4 */ |
---|
| 1141 | |
---|
| 1142 | |
---|
| 1143 | |
---|
| 1144 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1145 | void DecompTemplate( KPP_REAL A[], int Pivot[], int* ising ) |
---|
| 1146 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1147 | Template for the LU decomposition |
---|
| 1148 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1149 | { |
---|
| 1150 | *ising = KppDecomp ( A ); |
---|
| 1151 | /*~~~> Note: for a full matrix use Lapack: |
---|
| 1152 | DGETRF( KPP_NVAR, KPP_NVAR, A, KPP_NVAR, Pivot, ising ) */ |
---|
| 1153 | |
---|
| 1154 | Ndec++; |
---|
| 1155 | |
---|
| 1156 | } /* DecompTemplate */ |
---|
| 1157 | |
---|
| 1158 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1159 | void SolveTemplate( KPP_REAL A[], int Pivot[], KPP_REAL b[] ) |
---|
| 1160 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1161 | Template for the forward/backward substitution (using pre-computed LU decomposition) |
---|
| 1162 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1163 | { |
---|
| 1164 | KppSolve( A, b ); |
---|
| 1165 | /*~~~> Note: for a full matrix use Lapack: |
---|
| 1166 | NRHS = 1 |
---|
| 1167 | DGETRS( 'N', KPP_NVAR , NRHS, A, KPP_NVAR, Pivot, b, KPP_NVAR, INFO ) */ |
---|
| 1168 | |
---|
| 1169 | Nsol++; |
---|
| 1170 | |
---|
| 1171 | } /* SolveTemplate */ |
---|
| 1172 | |
---|
| 1173 | |
---|
| 1174 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1175 | void FunTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Ydot[] ) |
---|
| 1176 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1177 | Template for the ODE function call. |
---|
| 1178 | Updates the rate coefficients (and possibly the fixed species) at each call |
---|
| 1179 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1180 | { |
---|
| 1181 | KPP_REAL Told; |
---|
| 1182 | |
---|
| 1183 | Told = TIME; |
---|
| 1184 | TIME = T; |
---|
| 1185 | Update_SUN(); |
---|
| 1186 | Update_RCONST(); |
---|
| 1187 | Fun( Y, FIX, RCONST, Ydot ); |
---|
| 1188 | TIME = Told; |
---|
| 1189 | |
---|
| 1190 | Nfun++; |
---|
| 1191 | |
---|
| 1192 | } /* FunTemplate */ |
---|
| 1193 | |
---|
| 1194 | |
---|
| 1195 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1196 | void JacTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Jcb[] ) |
---|
| 1197 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1198 | Template for the ODE Jacobian call. |
---|
| 1199 | Updates the rate coefficients (and possibly the fixed species) at each call |
---|
| 1200 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1201 | { |
---|
| 1202 | /*~~~> Local variables */ |
---|
| 1203 | KPP_REAL Told; |
---|
| 1204 | |
---|
| 1205 | Told = TIME; |
---|
| 1206 | TIME = T ; |
---|
| 1207 | Update_SUN(); |
---|
| 1208 | Update_RCONST(); |
---|
| 1209 | Jac_SP( Y, FIX, RCONST, Jcb ); |
---|
| 1210 | TIME = Told; |
---|
| 1211 | |
---|
| 1212 | Njac++; |
---|
| 1213 | |
---|
| 1214 | } /* JacTemplate */ |
---|
| 1215 | |
---|
| 1216 | |
---|