[2696] | 1 | void INTEGRATE( double DT ) |
---|
| 2 | { |
---|
| 3 | KPP_REAL P_VAR[NVAR], D_VAR[NVAR], V1[NVAR], V2[NVAR]; |
---|
| 4 | int IsReject; |
---|
| 5 | KPP_REAL T, Tnext, STEP, STEPold, Told, SUP; |
---|
| 6 | KPP_REAL ERR, ERRold, ratio, factor, facmax, tmp; |
---|
| 7 | int i; |
---|
| 8 | |
---|
| 9 | T = TIME; |
---|
| 10 | Tnext = TIME + DT; |
---|
| 11 | STEP = STEPMIN; |
---|
| 12 | Told = T; |
---|
| 13 | SUP = 1e-14; |
---|
| 14 | IsReject = 0; |
---|
| 15 | ERR = .5; |
---|
| 16 | |
---|
| 17 | /* -- BELOW THIS LIMIT USE TAYLOR INSTEAD OF EXP --- */ |
---|
| 18 | |
---|
| 19 | while ( T < Tnext ) { |
---|
| 20 | |
---|
| 21 | T = Told + STEP; |
---|
| 22 | if ( T > Tnext ) { |
---|
| 23 | STEP = Tnext - Told; |
---|
| 24 | T = Tnext; |
---|
| 25 | } |
---|
| 26 | |
---|
| 27 | FSPLIT_VAR ( VAR, P_VAR, D_VAR ); |
---|
| 28 | |
---|
| 29 | for( i = 0; i < NVAR; i++ ) { |
---|
| 30 | if ( fabs(D_VAR[i]) > SUP ) { |
---|
| 31 | ratio = P_VAR[i] / D_VAR[i]; |
---|
| 32 | tmp = (KPP_REAL)exp( (double)(-D_VAR[i] * STEP * 0.5) ); |
---|
| 33 | V1[i] = tmp * tmp * (VAR[i] - ratio) + ratio; |
---|
| 34 | V2[i] = tmp * (VAR[i] - ratio) + ratio; |
---|
| 35 | } else { |
---|
| 36 | tmp = D_VAR[i] * STEP * 0.5; |
---|
| 37 | V1[i] = VAR[i] + P_VAR[i] * STEP * ( 1 - tmp * |
---|
| 38 | ( 1 - 2.0 / 3.0 * tmp ) ); |
---|
| 39 | V2[i] = VAR[i] + P_VAR[i] * 0.5 * STEP * ( 1 - 0.5 * tmp * |
---|
| 40 | ( 1 - 1.0 / 3.0 * tmp ) ); |
---|
| 41 | } |
---|
| 42 | } |
---|
| 43 | |
---|
| 44 | FSPLIT_VAR( V2, P_VAR, D_VAR ); |
---|
| 45 | |
---|
| 46 | for( i = 0; i < NVAR; i++ ) { |
---|
| 47 | if ( fabs(D_VAR[i]) > SUP ) { |
---|
| 48 | ratio = P_VAR[i] / D_VAR[i]; |
---|
| 49 | tmp = (KPP_REAL)exp( (double)(-D_VAR[i] * STEP * 0.5) ); |
---|
| 50 | V2[i] = tmp * (V2[i] - ratio) + ratio; |
---|
| 51 | } else { |
---|
| 52 | tmp = D_VAR[i] * STEP * 0.5; |
---|
| 53 | V2[i] = V2[i] + P_VAR[i] * 0.5 * STEP * ( 1 - 0.5 * tmp * |
---|
| 54 | ( 1 - 1.0 / 3.0 * tmp ) ); |
---|
| 55 | } |
---|
| 56 | } |
---|
| 57 | /* ==== Extrapolation and error estimation ======== */ |
---|
| 58 | |
---|
| 59 | ERRold=ERR; |
---|
| 60 | ERR=0.; |
---|
| 61 | for( i = 0; i < NVAR; i++ ) { |
---|
| 62 | V1[i] = 2.*V2[i] - V1[i]; |
---|
| 63 | tmp = (V2[i] - V1[i]) / (ATOL[i] + RTOL[i]*V2[i]); |
---|
| 64 | ERR = ERR + tmp*tmp; |
---|
| 65 | } |
---|
| 66 | ERR = sqrt(ERR/NVAR); |
---|
| 67 | STEPold = STEP; |
---|
| 68 | |
---|
| 69 | /* ===== choosing the stepsize ==================== */ |
---|
| 70 | |
---|
| 71 | factor = 0.9 / pow(ERR,0.35) * pow(ERRold,0.2); |
---|
| 72 | facmax = IsReject ? 1.0 : 5.0; |
---|
| 73 | |
---|
| 74 | factor = max( 0.2, min(factor,facmax) ); |
---|
| 75 | STEP = min( STEPMAX, max(STEPMIN,factor*STEP) ); |
---|
| 76 | |
---|
| 77 | /*================================================= */ |
---|
| 78 | |
---|
| 79 | if ( (ERR > 1) && (STEPold > STEPMIN) ) { |
---|
| 80 | T = Told; |
---|
| 81 | IsReject = 1; |
---|
| 82 | } else { |
---|
| 83 | IsReject = 0; |
---|
| 84 | Told = T; |
---|
| 85 | for( i = 0; i < NVAR; i++ ) |
---|
| 86 | VAR[i] = max( V1[i], 0.0 ); |
---|
| 87 | TIME = Tnext; |
---|
| 88 | } |
---|
| 89 | } |
---|
| 90 | } |
---|