| 1 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! |
|---|
| 2 | ! RADAU5 - Runge-Kutta method based on Radau-2A quadrature ! |
|---|
| 3 | ! (2 stages, order 5) ! |
|---|
| 4 | ! By default the code employs the KPP sparse linear algebra routines ! |
|---|
| 5 | ! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) ! |
|---|
| 6 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! |
|---|
| 7 | |
|---|
| 8 | MODULE KPP_ROOT_Integrator |
|---|
| 9 | |
|---|
| 10 | USE KPP_ROOT_Precision, ONLY: dp |
|---|
| 11 | USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO |
|---|
| 12 | USE KPP_ROOT_Jacobian, ONLY: LU_DIAG |
|---|
| 13 | USE KPP_ROOT_LinearAlgebra |
|---|
| 14 | |
|---|
| 15 | IMPLICIT NONE |
|---|
| 16 | PUBLIC |
|---|
| 17 | SAVE |
|---|
| 18 | |
|---|
| 19 | ! Statistics |
|---|
| 20 | INTEGER :: Nfun, Njac, Nstp, Nacc, Nrej, Ndec, Nsol, Nsng |
|---|
| 21 | |
|---|
| 22 | ! Method parameters |
|---|
| 23 | KPP_REAL :: Transf(3,3), TransfInv(3,3), & |
|---|
| 24 | rkA(3,3), rkB(3), rkC(3), rkE(3), & |
|---|
| 25 | rkGamma, rkAlpha, rkBeta |
|---|
| 26 | |
|---|
| 27 | ! description of the error numbers IERR |
|---|
| 28 | CHARACTER(LEN=50), PARAMETER, DIMENSION(-11:1) :: IERR_NAMES = (/ & |
|---|
| 29 | 'Matrix is repeatedly singular ', & ! -11 |
|---|
| 30 | 'Step size too small: T + 10*H = T or H < Roundoff ', & ! -10 |
|---|
| 31 | 'No of steps exceeds maximum bound ', & ! -9 |
|---|
| 32 | 'Tolerances are too small ', & ! -8 |
|---|
| 33 | 'Improper values for Qmin, Qmax ', & ! -7 |
|---|
| 34 | 'Newton stopping tolerance too small ', & ! -6 |
|---|
| 35 | 'Improper value for ThetaMin ', & ! -5 |
|---|
| 36 | 'Improper values for FacMin/FacMax/FacSafe/FacRej ', & ! -4 |
|---|
| 37 | 'Hmin/Hmax/Hstart must be positive ', & ! -3 |
|---|
| 38 | 'Improper value for maximal no of Newton iterations', & ! -2 |
|---|
| 39 | 'Improper value for maximal no of steps ', & ! -1 |
|---|
| 40 | ' ', & ! 0 (not used) |
|---|
| 41 | 'Success ' /) ! 1 |
|---|
| 42 | |
|---|
| 43 | CONTAINS |
|---|
| 44 | |
|---|
| 45 | ! ************************************************************************** |
|---|
| 46 | |
|---|
| 47 | SUBROUTINE INTEGRATE( TIN, TOUT, & |
|---|
| 48 | ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U ) |
|---|
| 49 | |
|---|
| 50 | USE KPP_ROOT_Parameters, ONLY: NVAR |
|---|
| 51 | USE KPP_ROOT_Global, ONLY: ATOL,RTOL,VAR |
|---|
| 52 | |
|---|
| 53 | IMPLICIT NONE |
|---|
| 54 | |
|---|
| 55 | KPP_REAL :: TIN ! TIN - Start Time |
|---|
| 56 | KPP_REAL :: TOUT ! TOUT - End Time |
|---|
| 57 | INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20) |
|---|
| 58 | KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20) |
|---|
| 59 | INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20) |
|---|
| 60 | KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20) |
|---|
| 61 | INTEGER, INTENT(OUT), OPTIONAL :: IERR_U |
|---|
| 62 | |
|---|
| 63 | KPP_REAL, SAVE :: H |
|---|
| 64 | INTEGER :: IERR |
|---|
| 65 | |
|---|
| 66 | KPP_REAL :: RCNTRL(20), RSTATUS(20) |
|---|
| 67 | INTEGER :: ICNTRL(20), ISTATUS(20) |
|---|
| 68 | INTEGER, SAVE :: Ntotal = 0 |
|---|
| 69 | |
|---|
| 70 | H =0.0_dp |
|---|
| 71 | |
|---|
| 72 | !~~~> fine-tune the integrator: |
|---|
| 73 | ICNTRL(:) = 0 |
|---|
| 74 | ICNTRL(2) = 0 ! 0=vector tolerances, 1=scalar tolerances |
|---|
| 75 | ICNTRL(5) = 8 ! Max no. of Newton iterations |
|---|
| 76 | ICNTRL(6) = 1 ! Starting values for Newton are interpolated (0) or zero (1) |
|---|
| 77 | ICNTRL(11) = 1 ! Gustaffson (1) or classic(2) controller |
|---|
| 78 | RCNTRL(1:20) = 0._dp |
|---|
| 79 | |
|---|
| 80 | !~~~> if optional parameters are given, and if they are >0, |
|---|
| 81 | ! then use them to overwrite default settings |
|---|
| 82 | IF (PRESENT(ICNTRL_U)) ICNTRL(1:20) = ICNTRL_U(1:20) |
|---|
| 83 | IF (PRESENT(RCNTRL_U)) RCNTRL(1:20) = RCNTRL_U(1:20) |
|---|
| 84 | |
|---|
| 85 | |
|---|
| 86 | CALL RADAU5( NVAR,TIN,TOUT,VAR,H, & |
|---|
| 87 | RTOL,ATOL, & |
|---|
| 88 | RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR ) |
|---|
| 89 | |
|---|
| 90 | !!$ Ntotal = Ntotal + Nstp |
|---|
| 91 | !!$ PRINT*,'NSTEPS=',Nstp,' (',Ntotal,')' |
|---|
| 92 | |
|---|
| 93 | Nfun = Nfun + ISTATUS(1) |
|---|
| 94 | Njac = Njac + ISTATUS(2) |
|---|
| 95 | Nstp = Nstp + ISTATUS(3) |
|---|
| 96 | Nacc = Nacc + ISTATUS(4) |
|---|
| 97 | Nrej = Nrej + ISTATUS(5) |
|---|
| 98 | Ndec = Ndec + ISTATUS(6) |
|---|
| 99 | Nsol = Nsol + ISTATUS(7) |
|---|
| 100 | Nsng = Nsng + ISTATUS(8) |
|---|
| 101 | |
|---|
| 102 | ! if optional parameters are given for output |
|---|
| 103 | ! use them to store information in them |
|---|
| 104 | IF (PRESENT(ISTATUS_U)) THEN |
|---|
| 105 | ISTATUS_U(:) = 0 |
|---|
| 106 | ISTATUS_U(1) = Nfun ! function calls |
|---|
| 107 | ISTATUS_U(2) = Njac ! jacobian calls |
|---|
| 108 | ISTATUS_U(3) = Nstp ! steps |
|---|
| 109 | ISTATUS_U(4) = Nacc ! accepted steps |
|---|
| 110 | ISTATUS_U(5) = Nrej ! rejected steps (except at the beginning) |
|---|
| 111 | ISTATUS_U(6) = Ndec ! LU decompositions |
|---|
| 112 | ISTATUS_U(7) = Nsol ! forward/backward substitutions |
|---|
| 113 | ENDIF |
|---|
| 114 | IF (PRESENT(RSTATUS_U)) THEN |
|---|
| 115 | RSTATUS_U(:) = 0. |
|---|
| 116 | RSTATUS_U(1) = TOUT ! final time |
|---|
| 117 | ENDIF |
|---|
| 118 | IF (PRESENT(IERR_U)) IERR_U = IERR |
|---|
| 119 | |
|---|
| 120 | ! mz_rs_20050716: IERR is returned to the user who then decides what to do |
|---|
| 121 | ! about it, i.e. either stop the run or ignore it. |
|---|
| 122 | !!$ IF (IERR < 0) THEN |
|---|
| 123 | !!$ PRINT *,'RADAU: Unsuccessful exit at T=', TIN,' (IERR=',IERR,')' |
|---|
| 124 | !!$ STOP |
|---|
| 125 | !!$ ENDIF |
|---|
| 126 | |
|---|
| 127 | END SUBROUTINE INTEGRATE |
|---|
| 128 | |
|---|
| 129 | SUBROUTINE RADAU5(N,T,Tend,Y,H,RelTol,AbsTol, & |
|---|
| 130 | RCNTRL,ICNTRL,RSTATUS,ISTATUS,IDID) |
|---|
| 131 | |
|---|
| 132 | !~~~>----------------------------------------------- |
|---|
| 133 | ! NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC) |
|---|
| 134 | ! SYSTEM OF FirstStep 0RDER ORDINARY DIFFERENTIAL EQUATIONS |
|---|
| 135 | ! M*Y'=F(T,Y). |
|---|
| 136 | ! THE SYSTEM CAN BE (LINEARLY) IMPLICIT (MASS-MATRIX M /= I) |
|---|
| 137 | ! OR EXPLICIT (M=I). |
|---|
| 138 | ! THE METHOD USED IS AN IMPLICIT RUNGE-KUTTA METHOD (RADAU IIA) |
|---|
| 139 | ! OF ORDER 5 WITH STEP SIZE CONTROL AND CONTINUOUS OUTPUT. |
|---|
| 140 | ! C.F. SECTION IV.8 |
|---|
| 141 | ! |
|---|
| 142 | ! AUTHORS: E. HAIRER AND G. WANNER |
|---|
| 143 | ! UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES |
|---|
| 144 | ! CH-1211 GENEVE 24, SWITZERLAND |
|---|
| 145 | ! E-MAIL: HAIRER@DIVSUN.UNIGE.CH, WANNER@DIVSUN.UNIGE.CH |
|---|
| 146 | ! |
|---|
| 147 | ! THIS CODE IS PART OF THE BOOK: |
|---|
| 148 | ! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL |
|---|
| 149 | ! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. |
|---|
| 150 | ! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14, |
|---|
| 151 | ! SPRINGER-VERLAG (1991) |
|---|
| 152 | ! |
|---|
| 153 | ! VERSION OF SEPTEMBER 30, 1995 |
|---|
| 154 | ! |
|---|
| 155 | ! INPUT PARAMETERS |
|---|
| 156 | ! ---------------- |
|---|
| 157 | ! N DIMENSION OF THE SYSTEM |
|---|
| 158 | ! |
|---|
| 159 | ! FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE |
|---|
| 160 | ! VALUE OF F(T,Y): |
|---|
| 161 | ! SUBROUTINE FCN(N,T,Y,F) |
|---|
| 162 | ! KPP_REAL T,Y(N),F(N) |
|---|
| 163 | ! F(1)=... ETC. |
|---|
| 164 | ! RPAR, IPAR (SEE BELOW) |
|---|
| 165 | ! |
|---|
| 166 | ! T INITIAL TIME VALUE |
|---|
| 167 | ! |
|---|
| 168 | ! Tend FINAL T-VALUE (Tend-T MAY BE POSITIVE OR NEGATIVE) |
|---|
| 169 | ! |
|---|
| 170 | ! Y(N) INITIAL VALUES FOR Y |
|---|
| 171 | ! |
|---|
| 172 | ! H INITIALL STEP SIZE GUESS; |
|---|
| 173 | ! FOR STIFF EQUATIONS WITH INITIALL TRANSIENT, |
|---|
| 174 | ! H=1.D0/(NORM OF F'), USUALLY 1.D-3 OR 1.D-5, IS GOOD. |
|---|
| 175 | ! THIS CHOICE IS NOT VERY IMPORTANT, THE STEP SIZE IS |
|---|
| 176 | ! QUICKLY ADAPTED. (IF H=0.D0, THE CODE PUTS H=1.D-6). |
|---|
| 177 | ! |
|---|
| 178 | ! RelTol,AbsTol RELATIVE AND ABSOLUTE ERROR TOLERANCES. |
|---|
| 179 | ! for ICNTRL(2) = 0: AbsTol, RelTol are N-dimensional vectors |
|---|
| 180 | ! = 1: AbsTol, RelTol are scalars |
|---|
| 181 | ! |
|---|
| 182 | ! ----- CONTINUOUS OUTPUT: ----- |
|---|
| 183 | ! DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION |
|---|
| 184 | ! FOR THE INTERVAL [Told,T] IS AVAILABLE THROUGH |
|---|
| 185 | ! THE FUNCTION |
|---|
| 186 | ! >>> CONTR5(I,S,CONT,LRC) <<< |
|---|
| 187 | ! WHICH PROVIDES AN APPROXIMATION TO THE I-TH |
|---|
| 188 | ! COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE |
|---|
| 189 | ! S SHOULD LIE IN THE INTERVAL [Told,T]. |
|---|
| 190 | ! DO NOT CHANGE THE ENTRIES OF CONT(LRC), IF THE |
|---|
| 191 | ! DENSE OUTPUT FUNCTION IS USED. |
|---|
| 192 | !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 193 | ! |
|---|
| 194 | !~~~> INPUT PARAMETERS: |
|---|
| 195 | ! |
|---|
| 196 | ! Note: For input parameters equal to zero the default values of the |
|---|
| 197 | ! corresponding variables are used. |
|---|
| 198 | ! |
|---|
| 199 | !~~~> Integer input parameters: |
|---|
| 200 | ! |
|---|
| 201 | ! ICNTRL(1) = not used |
|---|
| 202 | ! |
|---|
| 203 | ! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors |
|---|
| 204 | ! = 1: AbsTol, RelTol are scalars |
|---|
| 205 | ! |
|---|
| 206 | ! ICNTRL(3) = not used |
|---|
| 207 | ! |
|---|
| 208 | ! ICNTRL(4) -> maximum number of integration steps |
|---|
| 209 | ! For ICNTRL(4)=0 the default value of 10000 is used |
|---|
| 210 | ! |
|---|
| 211 | ! ICNTRL(5) -> maximum number of Newton iterations |
|---|
| 212 | ! For ICNTRL(5)=0 the default value of 8 is used |
|---|
| 213 | ! |
|---|
| 214 | ! ICNTRL(6) -> starting values of Newton iterations: |
|---|
| 215 | ! ICNTRL(6)=0 : starting values are obtained from |
|---|
| 216 | ! the extrapolated collocation solution |
|---|
| 217 | ! (the default) |
|---|
| 218 | ! ICNTRL(6)=1 : starting values are zero |
|---|
| 219 | ! |
|---|
| 220 | ! ICNTRL(11) -> switch for step size strategy |
|---|
| 221 | ! ICNTRL(8) == 1: mod. predictive controller (Gustafsson) |
|---|
| 222 | ! ICNTRL(8) == 2: classical step size control |
|---|
| 223 | ! the default value (for iwork(8)=0) is iwork(8)=1. |
|---|
| 224 | ! the choice iwork(8) == 1 seems to produce safer results; |
|---|
| 225 | ! for simple problems, the choice iwork(8) == 2 produces |
|---|
| 226 | ! often slightly faster runs |
|---|
| 227 | ! ( currently unused ) |
|---|
| 228 | ! |
|---|
| 229 | !~~~> Real input parameters: |
|---|
| 230 | ! |
|---|
| 231 | ! RCNTRL(1) -> not used |
|---|
| 232 | ! |
|---|
| 233 | ! RCNTRL(2) -> Hmax, upper bound for the integration step size |
|---|
| 234 | ! |
|---|
| 235 | ! RCNTRL(3) -> not used |
|---|
| 236 | ! |
|---|
| 237 | ! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2) |
|---|
| 238 | ! |
|---|
| 239 | ! RCNTRL(5) -> FacMax, upper bound on step increase factor (default=6) |
|---|
| 240 | ! |
|---|
| 241 | ! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections |
|---|
| 242 | ! (default=0.1) |
|---|
| 243 | ! |
|---|
| 244 | ! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller |
|---|
| 245 | ! than the predicted value (default=0.9) |
|---|
| 246 | ! |
|---|
| 247 | ! RCNTRL(8) -> ThetaMin. If Newton convergence rate smaller |
|---|
| 248 | ! than ThetaMin the Jacobian is not recomputed; |
|---|
| 249 | ! (default=0.001) |
|---|
| 250 | ! |
|---|
| 251 | ! RCNTRL(9) -> NewtonTol, stopping criterion for Newton's method |
|---|
| 252 | ! (default=0.03) |
|---|
| 253 | ! |
|---|
| 254 | ! RCNTRL(10) -> Qmin |
|---|
| 255 | ! |
|---|
| 256 | ! RCNTRL(11) -> Qmax. If Qmin < Hnew/Hold < Qmax, then the |
|---|
| 257 | ! step size is kept constant and the LU factorization |
|---|
| 258 | ! reused (default Qmin=1, Qmax=1.2) |
|---|
| 259 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 260 | ! |
|---|
| 261 | ! |
|---|
| 262 | ! OUTPUT PARAMETERS |
|---|
| 263 | ! ----------------- |
|---|
| 264 | ! T T-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED |
|---|
| 265 | ! (AFTER SUCCESSFUL RETURN T=Tend). |
|---|
| 266 | ! |
|---|
| 267 | ! Y(N) NUMERICAL SOLUTION AT T |
|---|
| 268 | ! |
|---|
| 269 | ! H PREDICTED STEP SIZE OF THE LastStep ACCEPTED STEP |
|---|
| 270 | ! |
|---|
| 271 | ! IDID REPORTS ON SUCCESSFULNESS UPON RETURN: |
|---|
| 272 | ! |
|---|
| 273 | ! ISTATUS(1) = No. of function calls |
|---|
| 274 | ! ISTATUS(2) = No. of jacobian calls |
|---|
| 275 | ! ISTATUS(3) = No. of steps |
|---|
| 276 | ! ISTATUS(4) = No. of accepted steps |
|---|
| 277 | ! ISTATUS(5) = No. of rejected steps (except at the beginning) |
|---|
| 278 | ! ISTATUS(6) = No. of LU decompositions |
|---|
| 279 | ! ISTATUS(7) = No. of forward/backward substitutions |
|---|
| 280 | ! ISTATUS(8) = No. of singular matrix decompositions |
|---|
| 281 | ! |
|---|
| 282 | ! RSTATUS(1) -> Texit, the time corresponding to the |
|---|
| 283 | ! computed Y upon return |
|---|
| 284 | ! RSTATUS(2) -> Hexit, last accepted step before exit |
|---|
| 285 | ! For multiple restarts, use Hexit as Hstart |
|---|
| 286 | ! in the subsequent run |
|---|
| 287 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 288 | |
|---|
| 289 | IMPLICIT NONE |
|---|
| 290 | |
|---|
| 291 | INTEGER :: N |
|---|
| 292 | KPP_REAL :: Y(N),AbsTol(N),RelTol(N),RCNTRL(20),RSTATUS(20) |
|---|
| 293 | INTEGER :: ICNTRL(20), ISTATUS(20) |
|---|
| 294 | LOGICAL :: StartNewton, Gustafsson |
|---|
| 295 | INTEGER :: IDID, ITOL |
|---|
| 296 | KPP_REAL :: H,Tend,T |
|---|
| 297 | |
|---|
| 298 | !~~~> Control arguments |
|---|
| 299 | INTEGER :: Max_no_steps, NewtonMaxit |
|---|
| 300 | KPP_REAL :: Hstart,Hmin,Hmax,Qmin,Qmax |
|---|
| 301 | KPP_REAL :: Roundoff, ThetaMin,TolNewton |
|---|
| 302 | KPP_REAL :: FacSafe,FacMin,FacMax,FacRej |
|---|
| 303 | !~~~> Local variables |
|---|
| 304 | INTEGER :: i |
|---|
| 305 | KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 |
|---|
| 306 | |
|---|
| 307 | !~~~> variables from the former COMMON block /CONRA5/ |
|---|
| 308 | ! KPP_REAL :: Tsol, Hsol |
|---|
| 309 | |
|---|
| 310 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 311 | ! SETTING THE PARAMETERS |
|---|
| 312 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 313 | Nfun=0 |
|---|
| 314 | Njac=0 |
|---|
| 315 | Nstp=0 |
|---|
| 316 | Nacc=0 |
|---|
| 317 | Nrej=0 |
|---|
| 318 | Ndec=0 |
|---|
| 319 | Nsol=0 |
|---|
| 320 | IDID = 0 |
|---|
| 321 | |
|---|
| 322 | !~~~> ICNTRL(1) - autonomous system - not used |
|---|
| 323 | !~~~> ITOL: 1 for vector and 0 for scalar AbsTol/RelTol |
|---|
| 324 | IF (ICNTRL(2) == 0) THEN |
|---|
| 325 | ITOL = 1 |
|---|
| 326 | ELSE |
|---|
| 327 | ITOL = 0 |
|---|
| 328 | END IF |
|---|
| 329 | !~~~> ICNTRL(3) - method selection - not used |
|---|
| 330 | !~~~> Max_no_steps: the maximal number of time steps |
|---|
| 331 | IF (ICNTRL(4) == 0) THEN |
|---|
| 332 | Max_no_steps = 10000 |
|---|
| 333 | ELSE |
|---|
| 334 | Max_no_steps=ICNTRL(4) |
|---|
| 335 | IF (Max_no_steps <= 0) THEN |
|---|
| 336 | WRITE(6,*) 'ICNTRL(4)=',ICNTRL(4) |
|---|
| 337 | CALL RAD_ErrorMsg(-1,T,ZERO,IDID) |
|---|
| 338 | END IF |
|---|
| 339 | END IF |
|---|
| 340 | !~~~> NewtonMaxit MAXIMAL NUMBER OF NEWTON ITERATIONS |
|---|
| 341 | IF (ICNTRL(5) == 0) THEN |
|---|
| 342 | NewtonMaxit = 8 |
|---|
| 343 | ELSE |
|---|
| 344 | NewtonMaxit=ICNTRL(5) |
|---|
| 345 | IF (NewtonMaxit <= 0) THEN |
|---|
| 346 | WRITE(6,*) 'ICNTRL(5)=',ICNTRL(5) |
|---|
| 347 | CALL RAD_ErrorMsg(-2,T,ZERO,IDID) |
|---|
| 348 | END IF |
|---|
| 349 | END IF |
|---|
| 350 | !~~~> StartNewton: Use extrapolation for starting values of Newton iterations |
|---|
| 351 | IF (ICNTRL(6) == 0) THEN |
|---|
| 352 | StartNewton = .TRUE. |
|---|
| 353 | ELSE |
|---|
| 354 | StartNewton = .FALSE. |
|---|
| 355 | END IF |
|---|
| 356 | !~~~> Gustafsson: step size controller |
|---|
| 357 | IF(ICNTRL(11) == 0)THEN |
|---|
| 358 | Gustafsson=.TRUE. |
|---|
| 359 | ELSE |
|---|
| 360 | Gustafsson=.FALSE. |
|---|
| 361 | END IF |
|---|
| 362 | |
|---|
| 363 | |
|---|
| 364 | !~~~> Roundoff SMALLEST NUMBER SATISFYING 1.0d0+Roundoff>1.0d0 |
|---|
| 365 | Roundoff=WLAMCH('E'); |
|---|
| 366 | |
|---|
| 367 | !~~~> RCNTRL(1) = Hmin - not used |
|---|
| 368 | Hmin = ZERO |
|---|
| 369 | !~~~> Hmax = maximal step size |
|---|
| 370 | IF (RCNTRL(2) == ZERO) THEN |
|---|
| 371 | Hmax=Tend-T |
|---|
| 372 | ELSE |
|---|
| 373 | Hmax=MAX(ABS(RCNTRL(7)),ABS(Tend-T)) |
|---|
| 374 | END IF |
|---|
| 375 | !~~~> RCNTRL(3) = Hstart - not used |
|---|
| 376 | Hstart = ZERO |
|---|
| 377 | !~~~> FacMin: lower bound on step decrease factor |
|---|
| 378 | IF(RCNTRL(4) == ZERO)THEN |
|---|
| 379 | FacMin = 0.2d0 |
|---|
| 380 | ELSE |
|---|
| 381 | FacMin = RCNTRL(4) |
|---|
| 382 | END IF |
|---|
| 383 | !~~~> FacMax: upper bound on step increase factor |
|---|
| 384 | IF(RCNTRL(5) == ZERO)THEN |
|---|
| 385 | FacMax = 8.D0 |
|---|
| 386 | ELSE |
|---|
| 387 | FacMax = RCNTRL(5) |
|---|
| 388 | END IF |
|---|
| 389 | !~~~> FacRej: step decrease factor after 2 consecutive rejections |
|---|
| 390 | IF(RCNTRL(6) == ZERO)THEN |
|---|
| 391 | FacRej = 0.1d0 |
|---|
| 392 | ELSE |
|---|
| 393 | FacRej = RCNTRL(6) |
|---|
| 394 | END IF |
|---|
| 395 | !~~~> FacSafe: by which the new step is slightly smaller |
|---|
| 396 | ! than the predicted value |
|---|
| 397 | IF (RCNTRL(7) == ZERO) THEN |
|---|
| 398 | FacSafe=0.9d0 |
|---|
| 399 | ELSE |
|---|
| 400 | FacSafe=RCNTRL(7) |
|---|
| 401 | END IF |
|---|
| 402 | IF ( (FacMax < ONE) .OR. (FacMin > ONE) .OR. & |
|---|
| 403 | (FacSafe <= 0.001D0) .OR. (FacSafe >= 1.0d0) ) THEN |
|---|
| 404 | WRITE(6,*)'RCNTRL(4:7)=',RCNTRL(4:7) |
|---|
| 405 | CALL RAD_ErrorMsg(-4,T,ZERO,IDID) |
|---|
| 406 | END IF |
|---|
| 407 | |
|---|
| 408 | !~~~> ThetaMin: decides whether the Jacobian should be recomputed |
|---|
| 409 | IF (RCNTRL(8) == ZERO) THEN |
|---|
| 410 | ThetaMin = 1.0d-3 |
|---|
| 411 | ELSE |
|---|
| 412 | ThetaMin=RCNTRL(8) |
|---|
| 413 | IF (ThetaMin <= 0.0d0 .OR. ThetaMin >= 1.0d0) THEN |
|---|
| 414 | WRITE(6,*) 'RCNTRL(8)=', RCNTRL(8) |
|---|
| 415 | CALL RAD_ErrorMsg(-5,T,ZERO,IDID) |
|---|
| 416 | END IF |
|---|
| 417 | END IF |
|---|
| 418 | !~~~> TolNewton: stopping crierion for Newton's method |
|---|
| 419 | IF (RCNTRL(9) == ZERO) THEN |
|---|
| 420 | TolNewton = 3.0d-2 |
|---|
| 421 | ELSE |
|---|
| 422 | TolNewton = RCNTRL(9) |
|---|
| 423 | IF (TolNewton <= Roundoff) THEN |
|---|
| 424 | WRITE(6,*) 'RCNTRL(9)=',RCNTRL(9) |
|---|
| 425 | CALL RAD_ErrorMsg(-6,T,ZERO,IDID) |
|---|
| 426 | END IF |
|---|
| 427 | END IF |
|---|
| 428 | !~~~> Qmin AND Qmax: IF Qmin < Hnew/Hold < Qmax, STEP SIZE = CONST. |
|---|
| 429 | IF (RCNTRL(10) == ZERO) THEN |
|---|
| 430 | Qmin=1.D0 |
|---|
| 431 | ELSE |
|---|
| 432 | Qmin=RCNTRL(10) |
|---|
| 433 | END IF |
|---|
| 434 | IF (RCNTRL(11) == ZERO) THEN |
|---|
| 435 | Qmax=1.2D0 |
|---|
| 436 | ELSE |
|---|
| 437 | Qmax=RCNTRL(11) |
|---|
| 438 | END IF |
|---|
| 439 | IF (Qmin > ONE .OR. Qmax < ONE) THEN |
|---|
| 440 | WRITE(6,*) 'RCNTRL(10:11)=',Qmin,Qmax |
|---|
| 441 | CALL RAD_ErrorMsg(-7,T,ZERO,IDID) |
|---|
| 442 | END IF |
|---|
| 443 | !~~~> Check if tolerances are reasonable |
|---|
| 444 | IF (ITOL == 0) THEN |
|---|
| 445 | IF (AbsTol(1) <= ZERO.OR.RelTol(1) <= 10.d0*Roundoff) THEN |
|---|
| 446 | WRITE (6,*) 'AbsTol/RelTol=',AbsTol,RelTol |
|---|
| 447 | CALL RAD_ErrorMsg(-8,T,ZERO,IDID) |
|---|
| 448 | END IF |
|---|
| 449 | ELSE |
|---|
| 450 | DO i=1,N |
|---|
| 451 | IF (AbsTol(i) <= ZERO.OR.RelTol(i) <= 10.d0*Roundoff) THEN |
|---|
| 452 | WRITE (6,*) 'AbsTol/RelTol(',i,')=',AbsTol(i),RelTol(i) |
|---|
| 453 | CALL RAD_ErrorMsg(-8,T,ZERO,IDID) |
|---|
| 454 | END IF |
|---|
| 455 | END DO |
|---|
| 456 | END IF |
|---|
| 457 | |
|---|
| 458 | !~~~> WHEN A FAIL HAS OCCURED, RETURN |
|---|
| 459 | IF (IDID < 0) RETURN |
|---|
| 460 | |
|---|
| 461 | |
|---|
| 462 | !~~~> CALL TO CORE INTEGRATOR ------------ |
|---|
| 463 | CALL RAD_Integrator( N,T,Y,Tend,Hmax,H,RelTol,AbsTol,ITOL,IDID, & |
|---|
| 464 | Max_no_steps,Roundoff,FacSafe,ThetaMin,TolNewton,Qmin,Qmax, & |
|---|
| 465 | NewtonMaxit,StartNewton,Gustafsson,FacMin,FacMax,FacRej ) |
|---|
| 466 | |
|---|
| 467 | ISTATUS(1)=Nfun |
|---|
| 468 | ISTATUS(2)=Njac |
|---|
| 469 | ISTATUS(3)=Nstp |
|---|
| 470 | ISTATUS(4)=Nacc |
|---|
| 471 | ISTATUS(5)=Nrej |
|---|
| 472 | ISTATUS(6)=Ndec |
|---|
| 473 | ISTATUS(7)=Nsol |
|---|
| 474 | ISTATUS(8)=Nsng |
|---|
| 475 | |
|---|
| 476 | ! END SUBROUTINE RADAU5 |
|---|
| 477 | CONTAINS ! INTERNAL PROCEDURES TO RADAU5 |
|---|
| 478 | |
|---|
| 479 | |
|---|
| 480 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 481 | SUBROUTINE RAD_Integrator( N,T,Y,Tend,Hmax,H,RelTol,AbsTol,ITOL,IDID, & |
|---|
| 482 | Max_no_steps,Roundoff,FacSafe,ThetaMin,TolNewton,Qmin,Qmax, & |
|---|
| 483 | NewtonMaxit,StartNewton,Gustafsson,FacMin,FacMax,FacRej ) |
|---|
| 484 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 485 | ! CORE INTEGRATOR FOR RADAU5 |
|---|
| 486 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 487 | |
|---|
| 488 | IMPLICIT NONE |
|---|
| 489 | INTEGER :: N |
|---|
| 490 | KPP_REAL Y(NVAR),Z1(NVAR),Z2(NVAR),Z3(NVAR),Y0(NVAR),& |
|---|
| 491 | SCAL(NVAR),F1(NVAR),F2(NVAR),F3(NVAR), & |
|---|
| 492 | CONT(N,4),AbsTol(NVAR),RelTol(NVAR) |
|---|
| 493 | |
|---|
| 494 | #ifdef FULL_ALGEBRA |
|---|
| 495 | KPP_REAL :: FJAC(NVAR,NVAR), E1(NVAR,NVAR) |
|---|
| 496 | DOUBLE COMPLEX :: E2(NVAR,NVAR) |
|---|
| 497 | #else |
|---|
| 498 | KPP_REAL :: FJAC(LU_NONZERO), E1(LU_NONZERO) |
|---|
| 499 | DOUBLE COMPLEX :: E2(LU_NONZERO) |
|---|
| 500 | #endif |
|---|
| 501 | |
|---|
| 502 | !~~~> Local variables |
|---|
| 503 | KPP_REAL :: TMP(NVAR), T, Tend, Tdirection, & |
|---|
| 504 | H, Hmax, HmaxN, Hacc, Hnew, Hopt, Hold, & |
|---|
| 505 | Fac, FacMin, Facmax, FacSafe, FacRej, FacGus, FacConv, & |
|---|
| 506 | Theta, ThetaMin, TolNewton, ERR, ERRACC, & |
|---|
| 507 | Qmin, Qmax, DYNO, Roundoff, & |
|---|
| 508 | AK, AK1, AK2, AK3, C3Q, & |
|---|
| 509 | Qnewton, DYTH, THQ, THQOLD, DYNOLD, & |
|---|
| 510 | DENOM, C1Q, C2Q, ALPHA, BETA, GAMMA, CFAC, ACONT3, QT |
|---|
| 511 | INTEGER :: IP1(NVAR),IP2(NVAR), ITOL, IDID, Max_no_steps, & |
|---|
| 512 | NewtonIter, NewtonMaxit, ISING |
|---|
| 513 | LOGICAL :: REJECT, FirstStep, FreshJac, LastStep, & |
|---|
| 514 | Gustafsson, StartNewton, NewtonDone |
|---|
| 515 | ! KPP_REAL, PARAMETER :: ONE = 1.0d0, ZERO = 0.0d0 |
|---|
| 516 | |
|---|
| 517 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 518 | ! INITIALISATIONS |
|---|
| 519 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 520 | |
|---|
| 521 | CALL RAD_Coefficients |
|---|
| 522 | |
|---|
| 523 | Tdirection=SIGN(1.D0,Tend-T) |
|---|
| 524 | HmaxN=MIN(ABS(Hmax),ABS(Tend-T)) |
|---|
| 525 | H=MIN(ABS(Hmin),ABS(Hstart)) |
|---|
| 526 | H=MIN(ABS(H),HmaxN) |
|---|
| 527 | IF (ABS(H) <= 10.D0*Roundoff) H=1.0D-6 |
|---|
| 528 | H=SIGN(H,Tdirection) |
|---|
| 529 | Hold=H |
|---|
| 530 | REJECT=.FALSE. |
|---|
| 531 | FirstStep=.TRUE. |
|---|
| 532 | LastStep=.FALSE. |
|---|
| 533 | FreshJac=.FALSE.; Theta=1.0d0 |
|---|
| 534 | IF ((T+H*1.0001D0-Tend)*Tdirection >= 0.D0) THEN |
|---|
| 535 | H=Tend-T |
|---|
| 536 | LastStep=.TRUE. |
|---|
| 537 | END IF |
|---|
| 538 | FacConv=1.D0 |
|---|
| 539 | CFAC=FacSafe*(1+2*NewtonMaxit) |
|---|
| 540 | Nsng=0 |
|---|
| 541 | ! Told=T |
|---|
| 542 | CALL RAD_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL) |
|---|
| 543 | CALL FUN_CHEM(T,Y,Y0) |
|---|
| 544 | |
|---|
| 545 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 546 | !~~~> Time loop begins |
|---|
| 547 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 548 | Tloop: DO WHILE ( (Tend-T)*Tdirection - Roundoff > ZERO ) |
|---|
| 549 | |
|---|
| 550 | !~~~> COMPUTE JACOBIAN MATRIX ANALYTICALLY |
|---|
| 551 | IF ( (.NOT.FreshJac) .AND. (Theta > ThetaMin) ) THEN |
|---|
| 552 | CALL JAC_CHEM(T,Y,FJAC) |
|---|
| 553 | FreshJac=.TRUE. |
|---|
| 554 | END IF |
|---|
| 555 | |
|---|
| 556 | !~~~> Compute the matrices E1 and E2 and their decompositions |
|---|
| 557 | GAMMA = rkGamma/H |
|---|
| 558 | ALPHA = rkAlpha/H |
|---|
| 559 | BETA = rkBeta/H |
|---|
| 560 | CALL RAD_DecompReal(N,FJAC,GAMMA,E1,IP1,ISING) |
|---|
| 561 | IF (ISING /= 0) THEN |
|---|
| 562 | Nsng=Nsng+1 |
|---|
| 563 | IF (Nsng >= 5) THEN |
|---|
| 564 | CALL RAD_ErrorMsg(-12,T,H,IDID); RETURN |
|---|
| 565 | END IF |
|---|
| 566 | H=H*0.5D0; REJECT=.TRUE.; LastStep=.FALSE. |
|---|
| 567 | CYCLE Tloop |
|---|
| 568 | END IF |
|---|
| 569 | CALL RAD_DecompCmplx(N,FJAC,ALPHA,BETA,E2,IP2,ISING) |
|---|
| 570 | IF (ISING /= 0) THEN |
|---|
| 571 | Nsng=Nsng+1 |
|---|
| 572 | IF (Nsng >= 5) THEN |
|---|
| 573 | CALL RAD_ErrorMsg(-12,T,H,IDID); RETURN |
|---|
| 574 | END IF |
|---|
| 575 | H=H*0.5D0; REJECT=.TRUE.; LastStep=.FALSE. |
|---|
| 576 | CYCLE Tloop |
|---|
| 577 | END IF |
|---|
| 578 | |
|---|
| 579 | 30 CONTINUE |
|---|
| 580 | Nstp=Nstp+1 |
|---|
| 581 | IF (Nstp > Max_no_steps) THEN |
|---|
| 582 | PRINT*,'Max number of time steps is ',Max_no_steps |
|---|
| 583 | CALL RAD_ErrorMsg(-9,T,H,IDID); RETURN |
|---|
| 584 | END IF |
|---|
| 585 | IF (0.1D0*ABS(H) <= ABS(T)*Roundoff) THEN |
|---|
| 586 | CALL RAD_ErrorMsg(-10,T,H,IDID); RETURN |
|---|
| 587 | END IF |
|---|
| 588 | |
|---|
| 589 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 590 | ! STARTING VALUES FOR NEWTON ITERATION |
|---|
| 591 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 592 | IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN |
|---|
| 593 | CALL Set2zero(N,Z1) |
|---|
| 594 | CALL Set2zero(N,Z2) |
|---|
| 595 | CALL Set2zero(N,Z3) |
|---|
| 596 | CALL Set2zero(N,F1) |
|---|
| 597 | CALL Set2zero(N,F2) |
|---|
| 598 | CALL Set2zero(N,F3) |
|---|
| 599 | ELSE |
|---|
| 600 | C3Q=H/Hold |
|---|
| 601 | C1Q=rkC(1)*C3Q |
|---|
| 602 | C2Q=rkC(2)*C3Q |
|---|
| 603 | DO i=1,N |
|---|
| 604 | AK1=CONT(i,2) |
|---|
| 605 | AK2=CONT(i,3) |
|---|
| 606 | AK3=CONT(i,4) |
|---|
| 607 | Z1(i)=C1Q*(AK1+(C1Q-rkC(2)+ONE)*(AK2+(C1Q-rkC(1)+ONE)*AK3)) |
|---|
| 608 | Z2(i)=C2Q*(AK1+(C2Q-rkC(2)+ONE)*(AK2+(C2Q-rkC(1)+ONE)*AK3)) |
|---|
| 609 | Z3(i)=C3Q*(AK1+(C3Q-rkC(2)+ONE)*(AK2+(C3Q-rkC(1)+ONE)*AK3)) |
|---|
| 610 | END DO |
|---|
| 611 | ! F(1,2,3) = TransfInv x Z(1,2,3) |
|---|
| 612 | CALL RAD_Transform(N,TransfInv,Z1,Z2,Z3,F1,F2,F3) |
|---|
| 613 | END IF |
|---|
| 614 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 615 | ! LOOP FOR THE SIMPLIFIED NEWTON ITERATION |
|---|
| 616 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 617 | |
|---|
| 618 | FacConv = MAX(FacConv,Roundoff)**0.8D0 |
|---|
| 619 | Theta=ABS(ThetaMin) |
|---|
| 620 | |
|---|
| 621 | NewtonLoop:DO NewtonIter = 1, NewtonMaxit |
|---|
| 622 | |
|---|
| 623 | !~~~> The right-hand side |
|---|
| 624 | DO i=1,N |
|---|
| 625 | TMP(i)=Y(i)+Z1(i) |
|---|
| 626 | END DO |
|---|
| 627 | CALL FUN_CHEM(T+rkC(1)*H,TMP,Z1) |
|---|
| 628 | DO i=1,N |
|---|
| 629 | TMP(i)=Y(i)+Z2(i) |
|---|
| 630 | END DO |
|---|
| 631 | CALL FUN_CHEM(T+rkC(2)*H,TMP,Z2) |
|---|
| 632 | DO i=1,N |
|---|
| 633 | TMP(i)=Y(i)+Z3(i) |
|---|
| 634 | END DO |
|---|
| 635 | CALL FUN_CHEM(T+rkC(3)*H,TMP,Z3) |
|---|
| 636 | |
|---|
| 637 | !~~~> Solve the linear systems |
|---|
| 638 | ! Z(1,2,3) = TransfInv x Z(1,2,3) |
|---|
| 639 | CALL RAD_Transform(N,TransfInv,Z1,Z2,Z3,Z1,Z2,Z3) |
|---|
| 640 | CALL RAD_Solve( N,FJAC,GAMMA,ALPHA,BETA,E1,E2, & |
|---|
| 641 | Z1,Z2,Z3,F1,F2,F3,CONT,IP1,IP2,ISING ) |
|---|
| 642 | Nsol=Nsol+1 |
|---|
| 643 | |
|---|
| 644 | DYNO=0.0d0 |
|---|
| 645 | DO i=1,N |
|---|
| 646 | DENOM=SCAL(i) |
|---|
| 647 | DYNO=DYNO+(Z1(i)/DENOM)**2+(Z2(i)/DENOM)**2+(Z3(i)/DENOM)**2 |
|---|
| 648 | END DO |
|---|
| 649 | DYNO=SQRT(DYNO/(3*N)) |
|---|
| 650 | |
|---|
| 651 | !~~~> Bad convergence or number of iterations too large |
|---|
| 652 | IF ( (NewtonIter > 1) .AND. (NewtonIter < NewtonMaxit) ) THEN |
|---|
| 653 | THQ=DYNO/DYNOLD |
|---|
| 654 | IF (NewtonIter == 2) THEN |
|---|
| 655 | Theta=THQ |
|---|
| 656 | ELSE |
|---|
| 657 | Theta=SQRT(THQ*THQOLD) |
|---|
| 658 | END IF |
|---|
| 659 | THQOLD=THQ |
|---|
| 660 | IF (Theta < 0.99d0) THEN |
|---|
| 661 | FacConv=Theta/(1.0d0-Theta) |
|---|
| 662 | DYTH=FacConv*DYNO*Theta**(NewtonMaxit-1-NewtonIter)/TolNewton |
|---|
| 663 | IF (DYTH >= 1.0d0) THEN |
|---|
| 664 | Qnewton=MAX(1.0D-4,MIN(20.0d0,DYTH)) |
|---|
| 665 | FAC=.8D0*Qnewton**(-1.0d0/(4.0d0+NewtonMaxit-1-NewtonIter)) |
|---|
| 666 | H=FAC*H |
|---|
| 667 | REJECT=.TRUE. |
|---|
| 668 | LastStep=.FALSE. |
|---|
| 669 | CYCLE Tloop |
|---|
| 670 | END IF |
|---|
| 671 | ELSE ! Non-convergence of Newton |
|---|
| 672 | H=H*0.5D0; REJECT=.TRUE.; LastStep=.FALSE. |
|---|
| 673 | CYCLE Tloop |
|---|
| 674 | END IF |
|---|
| 675 | END IF |
|---|
| 676 | DYNOLD=MAX(DYNO,Roundoff) |
|---|
| 677 | CALL WAXPY(N,ONE,Z1,1,F1,1) ! F1 <- F1 + Z1 |
|---|
| 678 | CALL WAXPY(N,ONE,Z2,1,F2,1) ! F2 <- F2 + Z2 |
|---|
| 679 | CALL WAXPY(N,ONE,Z3,1,F3,1) ! F3 <- F3 + Z3 |
|---|
| 680 | ! Z(1,2,3) = Transf x F(1,2,3) |
|---|
| 681 | CALL RAD_Transform(N,Transf,F1,F2,F3,Z1,Z2,Z3) |
|---|
| 682 | NewtonDone = (FacConv*DYNO <= TolNewton) |
|---|
| 683 | IF (NewtonDone) EXIT NewtonLoop |
|---|
| 684 | |
|---|
| 685 | END DO NewtonLoop |
|---|
| 686 | |
|---|
| 687 | IF (.NOT.NewtonDone) THEN |
|---|
| 688 | CALL RAD_ErrorMsg(-8,T,H,IDID); |
|---|
| 689 | H=H*0.5D0; REJECT=.TRUE.; LastStep=.FALSE. |
|---|
| 690 | CYCLE Tloop |
|---|
| 691 | END IF |
|---|
| 692 | |
|---|
| 693 | |
|---|
| 694 | !~~~> ERROR ESTIMATION |
|---|
| 695 | CALL RAD_ErrorEstimate(N,FJAC,H,Y0,Y,T, & |
|---|
| 696 | E1,Z1,Z2,Z3,IP1,SCAL,ERR, & |
|---|
| 697 | FirstStep,REJECT,GAMMA) |
|---|
| 698 | !~~~> COMPUTATION OF Hnew |
|---|
| 699 | Fac = ERR**(-0.25d0)* & |
|---|
| 700 | MIN(FacSafe,(NewtonIter+2*NewtonMaxit)/CFAC) |
|---|
| 701 | Fac = MIN(FacMax,MAX(FacMin,Fac)) |
|---|
| 702 | Hnew = Fac*H |
|---|
| 703 | |
|---|
| 704 | !~~~> IS THE ERROR SMALL ENOUGH ? |
|---|
| 705 | accept:IF (ERR < ONE) THEN !~~~> STEP IS ACCEPTED |
|---|
| 706 | FirstStep=.FALSE. |
|---|
| 707 | Nacc=Nacc+1 |
|---|
| 708 | IF (Gustafsson) THEN |
|---|
| 709 | !~~~> Predictive controller of Gustafsson |
|---|
| 710 | !~~~> Currently not implemented |
|---|
| 711 | IF (Nacc > 1) THEN |
|---|
| 712 | FacGus=FacSafe*(H/Hacc)*(ERR**2/ERRACC)**(-0.25d0) |
|---|
| 713 | FacGus=MIN(FacMax,MAX(FacMin,FacGus)) |
|---|
| 714 | Fac=MIN(Fac,FacGus) |
|---|
| 715 | Hnew=H*Fac |
|---|
| 716 | END IF |
|---|
| 717 | Hacc=H |
|---|
| 718 | ERRACC=MAX(1.0D-2,ERR) |
|---|
| 719 | END IF |
|---|
| 720 | ! Told = T |
|---|
| 721 | Hold = H |
|---|
| 722 | T=T+H |
|---|
| 723 | DO i=1,N |
|---|
| 724 | Y(i)=Y(i)+Z3(i) |
|---|
| 725 | CONT(i,2)=(Z2(i)-Z3(i))/(rkC(2)-ONE) |
|---|
| 726 | AK=(Z1(i)-Z2(i))/(rkC(1)-rkC(2)) |
|---|
| 727 | ACONT3=Z1(i)/rkC(1) |
|---|
| 728 | ACONT3=(AK-ACONT3)/rkC(2) |
|---|
| 729 | CONT(i,3)=(AK-CONT(i,2))/(rkC(1)-ONE) |
|---|
| 730 | CONT(i,4)=CONT(i,3)-ACONT3 |
|---|
| 731 | END DO |
|---|
| 732 | CALL RAD_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL) |
|---|
| 733 | FreshJac=.FALSE. |
|---|
| 734 | IF (LastStep) THEN |
|---|
| 735 | H=Hopt |
|---|
| 736 | IDID=1 |
|---|
| 737 | RETURN |
|---|
| 738 | END IF |
|---|
| 739 | CALL FUN_CHEM(T,Y,Y0) |
|---|
| 740 | Hnew=Tdirection*MIN(ABS(Hnew),HmaxN) |
|---|
| 741 | Hopt=Hnew |
|---|
| 742 | Hopt=MIN(H,Hnew) |
|---|
| 743 | IF (REJECT) Hnew=Tdirection*MIN(ABS(Hnew),ABS(H)) |
|---|
| 744 | REJECT=.FALSE. |
|---|
| 745 | IF ((T+Hnew/Qmin-Tend)*Tdirection >= 0.D0) THEN |
|---|
| 746 | H=Tend-T |
|---|
| 747 | LastStep=.TRUE. |
|---|
| 748 | ELSE |
|---|
| 749 | QT=Hnew/H |
|---|
| 750 | IF ( (Theta<=ThetaMin) .AND. (QT>=Qmin) & |
|---|
| 751 | .AND. (QT<=Qmax) ) GOTO 30 |
|---|
| 752 | H=Hnew |
|---|
| 753 | END IF |
|---|
| 754 | CYCLE Tloop |
|---|
| 755 | ELSE accept !~~~> STEP IS REJECTED |
|---|
| 756 | REJECT=.TRUE. |
|---|
| 757 | LastStep=.FALSE. |
|---|
| 758 | IF (FirstStep) THEN |
|---|
| 759 | H=H*FacRej |
|---|
| 760 | ELSE |
|---|
| 761 | H=Hnew |
|---|
| 762 | END IF |
|---|
| 763 | IF (Nacc >= 1) Nrej=Nrej+1 |
|---|
| 764 | CYCLE Tloop |
|---|
| 765 | END IF accept |
|---|
| 766 | |
|---|
| 767 | |
|---|
| 768 | END DO Tloop |
|---|
| 769 | |
|---|
| 770 | |
|---|
| 771 | END SUBROUTINE RAD_Integrator |
|---|
| 772 | |
|---|
| 773 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 774 | SUBROUTINE RAD_ErrorMsg(Code,T,H,IERR) |
|---|
| 775 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 776 | ! Handles all error messages |
|---|
| 777 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 778 | |
|---|
| 779 | IMPLICIT NONE |
|---|
| 780 | KPP_REAL, INTENT(IN) :: T, H |
|---|
| 781 | INTEGER, INTENT(IN) :: Code |
|---|
| 782 | INTEGER, INTENT(OUT) :: IERR |
|---|
| 783 | |
|---|
| 784 | IERR = Code |
|---|
| 785 | PRINT * , & |
|---|
| 786 | 'Forced exit from RADAU5 due to the following error:' |
|---|
| 787 | IF ((Code>=-11).AND.(Code<=-1)) THEN |
|---|
| 788 | PRINT *, IERR_NAMES(Code) |
|---|
| 789 | ELSE |
|---|
| 790 | PRINT *, 'Unknown Error code: ', Code |
|---|
| 791 | ENDIF |
|---|
| 792 | |
|---|
| 793 | PRINT *, "T=", T, "and H=", H |
|---|
| 794 | |
|---|
| 795 | END SUBROUTINE RAD_ErrorMsg |
|---|
| 796 | |
|---|
| 797 | |
|---|
| 798 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 799 | SUBROUTINE RAD_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL) |
|---|
| 800 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 801 | ! Handles all error messages |
|---|
| 802 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 803 | IMPLICIT NONE |
|---|
| 804 | INTEGER, INTENT(IN) :: N, ITOL |
|---|
| 805 | KPP_REAL, INTENT(IN) :: AbsTol(*), RelTol(*), Y(N) |
|---|
| 806 | KPP_REAL, INTENT(OUT) :: SCAL(N) |
|---|
| 807 | INTEGER :: i |
|---|
| 808 | |
|---|
| 809 | IF (ITOL==0) THEN |
|---|
| 810 | DO i=1,N |
|---|
| 811 | SCAL(i)=AbsTol(1)+RelTol(1)*ABS(Y(i)) |
|---|
| 812 | END DO |
|---|
| 813 | ELSE |
|---|
| 814 | DO i=1,N |
|---|
| 815 | SCAL(i)=AbsTol(i)+RelTol(i)*ABS(Y(i)) |
|---|
| 816 | END DO |
|---|
| 817 | END IF |
|---|
| 818 | |
|---|
| 819 | END SUBROUTINE RAD_ErrorScale |
|---|
| 820 | |
|---|
| 821 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 822 | SUBROUTINE RAD_Transform(N,Tr,Z1,Z2,Z3,F1,F2,F3) |
|---|
| 823 | !~~~> F = Tr x Z |
|---|
| 824 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 825 | IMPLICIT NONE |
|---|
| 826 | INTEGER :: N, i |
|---|
| 827 | KPP_REAL :: Tr(3,3),Z1(N),Z2(N),Z3(N),F1(N),F2(N),F3(N) |
|---|
| 828 | KPP_REAL :: x1, x2, x3 |
|---|
| 829 | DO i=1,N |
|---|
| 830 | x1 = Z1(i); x2 = Z2(i); x3 = Z3(i) |
|---|
| 831 | F1(i) = Tr(1,1)*x1 + Tr(1,2)*x2 + Tr(1,3)*x3 |
|---|
| 832 | F2(i) = Tr(2,1)*x1 + Tr(2,2)*x2 + Tr(2,3)*x3 |
|---|
| 833 | F3(i) = Tr(3,1)*x1 + Tr(3,2)*x2 + Tr(3,3)*x3 |
|---|
| 834 | END DO |
|---|
| 835 | END SUBROUTINE RAD_Transform |
|---|
| 836 | |
|---|
| 837 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 838 | SUBROUTINE RAD_Coefficients |
|---|
| 839 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 840 | IMPLICIT NONE |
|---|
| 841 | KPP_REAL :: s2,s3,s6,x1,x2,x3,x4,y1,y2,y3,y4,y5 |
|---|
| 842 | |
|---|
| 843 | s2 = SQRT(2.0d0); |
|---|
| 844 | s3 = SQRT(3.0d0); |
|---|
| 845 | s6 = SQRT(6.0d0); |
|---|
| 846 | x1 = 3.d0**(1.d0/3.d0); |
|---|
| 847 | x2 = 3.d0**(2.d0/3.d0); |
|---|
| 848 | x3 = 3.d0**(1.d0/6.d0); |
|---|
| 849 | x4 = 3.d0**(5.d0/6.d0); |
|---|
| 850 | |
|---|
| 851 | rkA(1,1) = 11.d0/45.d0-7.d0/360.d0*s6 |
|---|
| 852 | rkA(1,2) = 37.d0/225.d0-169.d0/1800.d0*s6 |
|---|
| 853 | rkA(1,3) = -2.d0/225.d0+s6/75 |
|---|
| 854 | rkA(2,1) = 37.d0/225.d0+169.d0/1800.d0*s6 |
|---|
| 855 | rkA(2,2) = 11.d0/45.d0+7.d0/360.d0*s6 |
|---|
| 856 | rkA(2,3) = -2.d0/225.d0-s6/75 |
|---|
| 857 | rkA(3,1) = 4.d0/9.d0-s6/36 |
|---|
| 858 | rkA(3,2) = 4.d0/9.d0+s6/36 |
|---|
| 859 | rkA(3,3) = 1.d0/9.d0 |
|---|
| 860 | |
|---|
| 861 | rkB(1) = 4.d0/9.d0-s6/36 |
|---|
| 862 | rkB(2) = 4.d0/9.d0+s6/36 |
|---|
| 863 | rkB(3) = 1.d0/9.d0 |
|---|
| 864 | |
|---|
| 865 | rkC(1) = 2.d0/5.d0-s6/10 |
|---|
| 866 | rkC(2) = 2.d0/5.d0+s6/10 |
|---|
| 867 | rkC(3) = 1.d0 |
|---|
| 868 | |
|---|
| 869 | ! Error estimation |
|---|
| 870 | rkE(1) = -(13.d0+7.d0*s6)/3.d0 |
|---|
| 871 | rkE(2) = (-13.d0+7.d0*s6)/3.d0 |
|---|
| 872 | rkE(3) = -1.d0/3.d0 |
|---|
| 873 | |
|---|
| 874 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 875 | !~~~> Diagonalize the RK matrix: |
|---|
| 876 | ! TransfInv * inv(rkA) * Transf = |
|---|
| 877 | ! | rkGamma 0 0 | |
|---|
| 878 | ! | 0 rkAlpha -rkBeta | |
|---|
| 879 | ! | 0 rkBeta rkAlpha | |
|---|
| 880 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 881 | rkGamma = 3-x1+x2 |
|---|
| 882 | rkAlpha = x1/2-x2/2+3 |
|---|
| 883 | rkBeta = -x4/2-3.d0/2.d0*x3 |
|---|
| 884 | |
|---|
| 885 | y1 = 36.d0/625.d0*s6 |
|---|
| 886 | y2 = 129.d0/2500.d0*x1 |
|---|
| 887 | y3 = 111.d0/2500.d0*x3*s2 |
|---|
| 888 | Transf(1,1) = -31.d0/1250.d0*s6*x1+37.d0/1250.d0*s6*x2-y1 & |
|---|
| 889 | +129.d0/1250.d0*x1-33.d0/1250.d0*x2+49.d0/625.d0 |
|---|
| 890 | Transf(1,2) = -y1-y2-y3 & |
|---|
| 891 | +31.d0/2500.d0*x4*s2+33.d0/2500.d0*x2+49.d0/625.d0 |
|---|
| 892 | Transf(1,3) = 3.d0/2500.d0*x3*(-33-43*x2+31*x3*s2+37*s3*s2) |
|---|
| 893 | Transf(2,1) = 31.d0/1250.d0*s6*x1-37.d0/1250.d0*s6*x2+y1 & |
|---|
| 894 | +129.d0/1250.d0*x1-33.d0/1250.d0*x2+49.d0/625.d0 |
|---|
| 895 | Transf(2,2) = y1-y2+y3& |
|---|
| 896 | -31.d0/2500.d0*x4*s2+33.d0/2500.d0*x2+49.d0/625.d0 |
|---|
| 897 | Transf(2,3) = -3.d0/2500.d0*x3*(33+43*x2+31*x3*s2+37*s3*s2) |
|---|
| 898 | Transf(3,1) = 1.d0 |
|---|
| 899 | Transf(3,2) = 1.d0 |
|---|
| 900 | Transf(3,3) = 0.d0 |
|---|
| 901 | |
|---|
| 902 | y1 = 11.d0/36.d0*x3*s2 + 43.d0/108.d0*x4*s2 |
|---|
| 903 | y2 = 11.d0/36.d0*s2*x2 - 43.d0/36.d0*s2*x1 |
|---|
| 904 | y3 = 31.d0/54.d0*x1 + 37.d0/54.d0*x2 |
|---|
| 905 | y4 = 31.d0/54.d0*x4-37.d0/18.d0*x3 |
|---|
| 906 | y5 = -x2/27+5.d0/27.d0*x1 |
|---|
| 907 | TransfInv(1,1) = y1 + y3 |
|---|
| 908 | TransfInv(1,2) = -y1 + y3 |
|---|
| 909 | TransfInv(1,3) = y5 + 1.d0/3.d0 |
|---|
| 910 | TransfInv(2,1) = -y1 - y3 |
|---|
| 911 | TransfInv(2,2) = y1 - y3 |
|---|
| 912 | TransfInv(2,3) = -y5 + 2.d0/3.d0 |
|---|
| 913 | TransfInv(3,1) = y4 - y2 |
|---|
| 914 | TransfInv(3,2) = y4 + y2 |
|---|
| 915 | TransfInv(3,3) = x3/9+5.d0/27.d0*x4 |
|---|
| 916 | |
|---|
| 917 | END SUBROUTINE RAD_Coefficients |
|---|
| 918 | |
|---|
| 919 | |
|---|
| 920 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 921 | SUBROUTINE RAD_DecompReal(N,FJAC,GAMMA,E1,IP1,ISING) |
|---|
| 922 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 923 | IMPLICIT NONE |
|---|
| 924 | |
|---|
| 925 | INTEGER :: N, ISING |
|---|
| 926 | KPP_REAL :: GAMMA |
|---|
| 927 | #ifdef FULL_ALGEBRA |
|---|
| 928 | KPP_REAL :: FJAC(NVAR,NVAR),E1(NVAR,NVAR) |
|---|
| 929 | #else |
|---|
| 930 | KPP_REAL :: FJAC(LU_NONZERO),E1(LU_NONZERO) |
|---|
| 931 | #endif |
|---|
| 932 | INTEGER :: IP1(N), i, j |
|---|
| 933 | |
|---|
| 934 | #ifdef FULL_ALGEBRA |
|---|
| 935 | DO j=1,N |
|---|
| 936 | DO i=1,N |
|---|
| 937 | E1(i,j)=-FJAC(i,j) |
|---|
| 938 | END DO |
|---|
| 939 | E1(j,j)=E1(j,j)+GAMMA |
|---|
| 940 | END DO |
|---|
| 941 | CALL DGETRF(N,N,E1,N,IP1,ISING) |
|---|
| 942 | #else |
|---|
| 943 | DO i=1,LU_NONZERO |
|---|
| 944 | E1(i)=-FJAC(i) |
|---|
| 945 | END DO |
|---|
| 946 | DO i=1,NVAR |
|---|
| 947 | j = LU_DIAG(i); E1(j)=E1(j)+GAMMA |
|---|
| 948 | END DO |
|---|
| 949 | CALL KppDecomp(E1,ISING) |
|---|
| 950 | #endif |
|---|
| 951 | |
|---|
| 952 | END SUBROUTINE RAD_DecompReal |
|---|
| 953 | |
|---|
| 954 | |
|---|
| 955 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 956 | SUBROUTINE RAD_DecompCmplx(N,FJAC,ALPHA,BETA,E2,IP2,ISING) |
|---|
| 957 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 958 | IMPLICIT NONE |
|---|
| 959 | INTEGER :: N, ISING |
|---|
| 960 | #ifdef FULL_ALGEBRA |
|---|
| 961 | KPP_REAL :: FJAC(N,N) |
|---|
| 962 | DOUBLE COMPLEX :: E2(N,N) |
|---|
| 963 | #else |
|---|
| 964 | KPP_REAL :: FJAC(LU_NONZERO) |
|---|
| 965 | DOUBLE COMPLEX :: E2(LU_NONZERO) |
|---|
| 966 | #endif |
|---|
| 967 | KPP_REAL :: ALPHA, BETA |
|---|
| 968 | INTEGER :: IP2(N), i, j |
|---|
| 969 | |
|---|
| 970 | #ifdef FULL_ALGEBRA |
|---|
| 971 | DO j=1,N |
|---|
| 972 | DO i=1,N |
|---|
| 973 | E2(i,j) = DCMPLX( -FJAC(i,j), 0.0d0 ) |
|---|
| 974 | END DO |
|---|
| 975 | E2(j,j) = E2(j,j) + DCMPLX( ALPHA, BETA ) |
|---|
| 976 | END DO |
|---|
| 977 | CALL ZGETRF(N,N,E2,N,IP2,ISING) |
|---|
| 978 | #else |
|---|
| 979 | DO i=1,LU_NONZERO |
|---|
| 980 | E2(i) = DCMPLX( -FJAC(i), 0.0d0 ) |
|---|
| 981 | END DO |
|---|
| 982 | DO i=1,NVAR |
|---|
| 983 | j = LU_DIAG(i); E2(j)=E2(j)+DCMPLX( ALPHA, BETA ) |
|---|
| 984 | END DO |
|---|
| 985 | CALL KppDecompCmplx(E2,ISING) |
|---|
| 986 | #endif |
|---|
| 987 | Ndec=Ndec+1 |
|---|
| 988 | |
|---|
| 989 | END SUBROUTINE RAD_DecompCmplx |
|---|
| 990 | |
|---|
| 991 | |
|---|
| 992 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 993 | SUBROUTINE RAD_Solve(N,FJAC,GAMMA,ALPHA,BETA,E1,E2,& |
|---|
| 994 | Z1,Z2,Z3,F1,F2,F3,CONT,IP1,IP2,ISING) |
|---|
| 995 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 996 | IMPLICIT NONE |
|---|
| 997 | INTEGER :: N,IP1(NVAR),IP2(NVAR),ISING |
|---|
| 998 | #ifdef FULL_ALGEBRA |
|---|
| 999 | KPP_REAL :: FJAC(NVAR,NVAR), E1(NVAR,NVAR) |
|---|
| 1000 | DOUBLE COMPLEX :: E2(NVAR,NVAR) |
|---|
| 1001 | #else |
|---|
| 1002 | KPP_REAL :: FJAC(LU_NONZERO), E1(LU_NONZERO) |
|---|
| 1003 | DOUBLE COMPLEX :: E2(LU_NONZERO) |
|---|
| 1004 | #endif |
|---|
| 1005 | KPP_REAL :: Z1(N),Z2(N),Z3(N), & |
|---|
| 1006 | F1(N),F2(N),F3(N),CONT(N), & |
|---|
| 1007 | GAMMA,ALPHA,BETA |
|---|
| 1008 | DOUBLE COMPLEX :: BC(N) |
|---|
| 1009 | INTEGER :: i,j |
|---|
| 1010 | KPP_REAL :: S2, S3 |
|---|
| 1011 | ! |
|---|
| 1012 | DO i=1,N |
|---|
| 1013 | S2=-F2(i) |
|---|
| 1014 | S3=-F3(i) |
|---|
| 1015 | Z1(i)=Z1(i)-F1(i)*GAMMA |
|---|
| 1016 | Z2(i)=Z2(i)+S2*ALPHA-S3*BETA |
|---|
| 1017 | Z3(i)=Z3(i)+S3*ALPHA+S2*BETA |
|---|
| 1018 | END DO |
|---|
| 1019 | #ifdef FULL_ALGEBRA |
|---|
| 1020 | CALL DGETRS ('N',N,1,E1,N,IP1,Z1,N,ISING) |
|---|
| 1021 | #else |
|---|
| 1022 | CALL KppSolve (E1,Z1) |
|---|
| 1023 | #endif |
|---|
| 1024 | |
|---|
| 1025 | DO j=1,N |
|---|
| 1026 | BC(j) = DCMPLX(Z2(j),Z3(j)) |
|---|
| 1027 | END DO |
|---|
| 1028 | #ifdef FULL_ALGEBRA |
|---|
| 1029 | CALL ZGETRS ('N',N,1,E2,N,IP2,BC,N,ISING) |
|---|
| 1030 | #else |
|---|
| 1031 | CALL KppSolveCmplx (E2,BC) |
|---|
| 1032 | #endif |
|---|
| 1033 | DO j=1,N |
|---|
| 1034 | Z2(j) = DBLE( BC(j) ) |
|---|
| 1035 | Z3(j) = AIMAG( BC(j) ) |
|---|
| 1036 | END DO |
|---|
| 1037 | |
|---|
| 1038 | END SUBROUTINE RAD_Solve |
|---|
| 1039 | |
|---|
| 1040 | |
|---|
| 1041 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1042 | SUBROUTINE RAD_ErrorEstimate(N,FJAC,H,Y0,Y,T,& |
|---|
| 1043 | E1,Z1,Z2,Z3,IP1,SCAL,ERR, & |
|---|
| 1044 | FirstStep,REJECT,GAMMA) |
|---|
| 1045 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1046 | IMPLICIT NONE |
|---|
| 1047 | |
|---|
| 1048 | INTEGER :: N |
|---|
| 1049 | #ifdef FULL_ALGEBRA |
|---|
| 1050 | KPP_REAL :: FJAC(NVAR,NVAR), E1(NVAR,NVAR) |
|---|
| 1051 | INTEGER :: ISING |
|---|
| 1052 | #else |
|---|
| 1053 | KPP_REAL :: FJAC(LU_NONZERO), E1(LU_NONZERO) |
|---|
| 1054 | #endif |
|---|
| 1055 | KPP_REAL :: SCAL(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N), & |
|---|
| 1056 | Y0(N),Y(N),TMP(N),T,H,GAMMA |
|---|
| 1057 | INTEGER :: IP1(N), i |
|---|
| 1058 | LOGICAL FirstStep,REJECT |
|---|
| 1059 | KPP_REAL :: HEE1,HEE2,HEE3,ERR |
|---|
| 1060 | |
|---|
| 1061 | HEE1 = rkE(1)/H |
|---|
| 1062 | HEE2 = rkE(2)/H |
|---|
| 1063 | HEE3 = rkE(3)/H |
|---|
| 1064 | |
|---|
| 1065 | DO i=1,N |
|---|
| 1066 | F2(i)=HEE1*Z1(i)+HEE2*Z2(i)+HEE3*Z3(i) |
|---|
| 1067 | TMP(i)=F2(i)+Y0(i) |
|---|
| 1068 | END DO |
|---|
| 1069 | |
|---|
| 1070 | #ifdef FULL_ALGEBRA |
|---|
| 1071 | CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,ISING) |
|---|
| 1072 | #else |
|---|
| 1073 | CALL KppSolve (E1, TMP) |
|---|
| 1074 | #endif |
|---|
| 1075 | |
|---|
| 1076 | ERR=0.D0 |
|---|
| 1077 | DO i=1,N |
|---|
| 1078 | ERR=ERR+(TMP(i)/SCAL(i))**2 |
|---|
| 1079 | END DO |
|---|
| 1080 | ERR=MAX(SQRT(ERR/N),1.D-10) |
|---|
| 1081 | ! |
|---|
| 1082 | IF (ERR < 1.D0) RETURN |
|---|
| 1083 | firej:IF (FirstStep.OR.REJECT) THEN |
|---|
| 1084 | DO i=1,N |
|---|
| 1085 | TMP(i)=Y(i)+TMP(i) |
|---|
| 1086 | END DO |
|---|
| 1087 | CALL FUN_CHEM(T,TMP,F1) |
|---|
| 1088 | DO i=1,N |
|---|
| 1089 | TMP(i)=F1(i)+F2(i) |
|---|
| 1090 | END DO |
|---|
| 1091 | |
|---|
| 1092 | #ifdef FULL_ALGEBRA |
|---|
| 1093 | CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,0) |
|---|
| 1094 | #else |
|---|
| 1095 | CALL KppSolve (E1, TMP) |
|---|
| 1096 | #endif |
|---|
| 1097 | ERR=0.D0 |
|---|
| 1098 | DO i=1,N |
|---|
| 1099 | ERR=ERR+(TMP(i)/SCAL(i))**2 |
|---|
| 1100 | END DO |
|---|
| 1101 | ERR=MAX(SQRT(ERR/N),1.0d-10) |
|---|
| 1102 | END IF firej |
|---|
| 1103 | |
|---|
| 1104 | END SUBROUTINE RAD_ErrorEstimate |
|---|
| 1105 | |
|---|
| 1106 | !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1107 | ! KPP_REAL FUNCTION CONTR5(I,N,T,CONT) |
|---|
| 1108 | !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1109 | !! THIS FUNCTION CAN BE USED FOR CONTINUOUS OUTPUT. IT PROVIDES AN |
|---|
| 1110 | !! APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT T. |
|---|
| 1111 | !! IT GIVES THE VALUE OF THE COLLOCATION POLYNOMIAL, DEFINED FOR |
|---|
| 1112 | !! THE STEP SUCCESSFULLY COMPUTED STEP (BY RADAU5). |
|---|
| 1113 | !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1114 | ! IMPLICIT NONE |
|---|
| 1115 | ! INTEGER :: I, N |
|---|
| 1116 | ! KPP_REAL :: T, CONT(N,4) |
|---|
| 1117 | ! KPP_REAL :: S |
|---|
| 1118 | ! KPP_REAL, PARAMETER :: ONE = 1.0d0 |
|---|
| 1119 | ! S=(T-Tsol)/Hsol |
|---|
| 1120 | ! CONTR5=CONT(i,1)+S* & |
|---|
| 1121 | ! (CONT(i,2)+(S-rkC(2)+ONE)*(CONT(i,3)+(S-rkC(1)+ONE)*CONT(i,4))) |
|---|
| 1122 | ! END FUNCTION CONTR5 |
|---|
| 1123 | |
|---|
| 1124 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1125 | END SUBROUTINE RADAU5 ! AND ALL ITS INTERNAL PROCEDURES |
|---|
| 1126 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1127 | |
|---|
| 1128 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1129 | SUBROUTINE FUN_CHEM(T, V, FCT) |
|---|
| 1130 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1131 | |
|---|
| 1132 | USE KPP_ROOT_Parameters |
|---|
| 1133 | USE KPP_ROOT_Global |
|---|
| 1134 | USE KPP_ROOT_Function, ONLY: Fun |
|---|
| 1135 | USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO |
|---|
| 1136 | |
|---|
| 1137 | IMPLICIT NONE |
|---|
| 1138 | |
|---|
| 1139 | KPP_REAL :: V(NVAR), FCT(NVAR) |
|---|
| 1140 | KPP_REAL :: T, Told |
|---|
| 1141 | |
|---|
| 1142 | !Told = TIME |
|---|
| 1143 | !TIME = T |
|---|
| 1144 | !CALL Update_SUN() |
|---|
| 1145 | !CALL Update_RCONST() |
|---|
| 1146 | !CALL Update_PHOTO() |
|---|
| 1147 | !TIME = Told |
|---|
| 1148 | |
|---|
| 1149 | CALL Fun(V, FIX, RCONST, FCT) |
|---|
| 1150 | |
|---|
| 1151 | Nfun=Nfun+1 |
|---|
| 1152 | |
|---|
| 1153 | END SUBROUTINE FUN_CHEM |
|---|
| 1154 | |
|---|
| 1155 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1156 | SUBROUTINE JAC_CHEM (T, V, JF) |
|---|
| 1157 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1158 | |
|---|
| 1159 | USE KPP_ROOT_Parameters |
|---|
| 1160 | USE KPP_ROOT_Global |
|---|
| 1161 | USE KPP_ROOT_JacobianSP |
|---|
| 1162 | USE KPP_ROOT_Jacobian, ONLY: Jac_SP |
|---|
| 1163 | USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO |
|---|
| 1164 | |
|---|
| 1165 | IMPLICIT NONE |
|---|
| 1166 | |
|---|
| 1167 | KPP_REAL :: V(NVAR), T, Told |
|---|
| 1168 | #ifdef FULL_ALGEBRA |
|---|
| 1169 | KPP_REAL :: JV(LU_NONZERO), JF(NVAR,NVAR) |
|---|
| 1170 | INTEGER :: i, j |
|---|
| 1171 | #else |
|---|
| 1172 | KPP_REAL :: JF(LU_NONZERO) |
|---|
| 1173 | #endif |
|---|
| 1174 | |
|---|
| 1175 | !Told = TIME |
|---|
| 1176 | !TIME = T |
|---|
| 1177 | !CALL Update_SUN() |
|---|
| 1178 | !CALL Update_RCONST() |
|---|
| 1179 | !CALL Update_PHOTO() |
|---|
| 1180 | !TIME = Told |
|---|
| 1181 | |
|---|
| 1182 | #ifdef FULL_ALGEBRA |
|---|
| 1183 | CALL Jac_SP(V, FIX, RCONST, JV) |
|---|
| 1184 | DO j=1,NVAR |
|---|
| 1185 | DO i=1,NVAR |
|---|
| 1186 | JF(i,j) = 0.0d0 |
|---|
| 1187 | END DO |
|---|
| 1188 | END DO |
|---|
| 1189 | DO i=1,LU_NONZERO |
|---|
| 1190 | JF(LU_IROW(i),LU_ICOL(i)) = JV(i) |
|---|
| 1191 | END DO |
|---|
| 1192 | #else |
|---|
| 1193 | CALL Jac_SP(V, FIX, RCONST, JF) |
|---|
| 1194 | #endif |
|---|
| 1195 | |
|---|
| 1196 | Njac=Njac+1 |
|---|
| 1197 | |
|---|
| 1198 | END SUBROUTINE JAC_CHEM |
|---|
| 1199 | |
|---|
| 1200 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1201 | |
|---|
| 1202 | END MODULE KPP_ROOT_Integrator |
|---|
