?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge remote-tracking branch 'origin/release-v4.6.1'
[WRF.git] / chem / KPP / kpp / kpp-2.1 / int / rosenbrock_adj.f90
blobfe9567fee89bff14879a2680e00e677029c1054e
1 MODULE KPP_ROOT_Integrator
2
3 USE KPP_ROOT_Precision
4 USE KPP_ROOT_Parameters
5 USE KPP_ROOT_Global
6 USE KPP_ROOT_LinearAlgebra
7 USE KPP_ROOT_Rates
8 USE KPP_ROOT_Function
9 USE KPP_ROOT_Jacobian
10 USE KPP_ROOT_Hessian
11 USE KPP_ROOT_Util
12
13 IMPLICIT NONE
14 PUBLIC
15 SAVE
16 !~~~> Statistics on the work performed by the Rosenbrock method
17 INTEGER :: Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
18 INTEGER, PARAMETER :: ifun=11, ijac=12, istp=13, iacc=14, &
19 irej=15, idec=16, isol=17, isng=18, &
20 itexit=11,ihexit=12
21 !~~~> Types of Adjoints Implemented
22 INTEGER, PARAMETER :: Adj_none = 1, Adj_discrete = 2, &
23 Adj_continuous = 3, Adj_simple_continuous = 4
24 KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
25 !~~~> Checkpoints in memory
26 INTEGER, PARAMETER :: bufsize = 1500
27 INTEGER :: stack_ptr = 0 ! last written entry
28 KPP_REAL, DIMENSION(:), POINTER :: buf_H, buf_T
29 KPP_REAL, DIMENSION(:,:), POINTER :: buf_Y, buf_K, buf_J
30 KPP_REAL, DIMENSION(:,:), POINTER :: buf_dY, buf_d2Y
31
32 CONTAINS ! Functions in the module KPP_ROOT_Integrator
33
34
35 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
36 SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, &
37 ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U )
38 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
39
40 IMPLICIT NONE
41
42 !~~~> Y - Concentrations
43 KPP_REAL :: Y(NVAR)
44 !~~~> NADJ - No. of cost functionals for which adjoints
45 ! are evaluated simultaneously
46 ! If single cost functional is considered (like in
47 ! most applications) simply set NADJ = 1
48 INTEGER NADJ
49 !~~~> Lambda - Sensitivities of concentrations
50 ! Note: Lambda (1:NVAR,j) contains sensitivities of
51 ! the j-th cost functional w.r.t. Y(1:NVAR), j=1...NADJ
52 KPP_REAL :: Lambda(NVAR,NADJ)
53 KPP_REAL, INTENT(IN) :: TIN ! TIN - Start Time
54 KPP_REAL, INTENT(IN) :: TOUT ! TOUT - End Time
55 !~~~> Optional input parameters and statistics
56 INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20)
57 KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20)
58 INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
59 KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
60
61 INTEGER, SAVE :: N_stp, N_acc, N_rej, N_sng, IERR
62 INTEGER :: i
63 KPP_REAL :: RCNTRL(20), RSTATUS(20)
64 INTEGER :: ICNTRL(20), ISTATUS(20)
65
66
67 ICNTRL(1:20) = 0
68 RCNTRL(1:20) = 0.0_dp
69 ISTATUS(1:20) = 0
70 RSTATUS(1:20) = 0.0_dp
71
72
73 ICNTRL(1) = 0 ! 0 = non-autonomous, 1 = autonomous
74 ICNTRL(2) = 1 ! 0 = scalar, 1 = vector tolerances
75 RCNTRL(3) = STEPMIN ! starting step
76 ICNTRL(4) = 5 ! choice of the method for forward integration
77 ICNTRL(5) = 2 ! 1=none, 2=discrete, 3=full continuous, 4=simplified continuous adjoint
78 ICNTRL(6) = 1 ! choice of the method for continuous adjoint
79
80 ! Tighter tolerances, especially atol, are needed for the full continuous adjoint
81 ! (Atol on sensitivities is different than on concentrations)
82 ! CADJ_ATOL(1:NVAR) = 1.0d-5
83 ! CADJ_RTOL(1:NVAR) = 1.0d-4
84
85 ! if optional parameters are given, and if they are >=0, then they overwrite default settings
86 IF (PRESENT(ICNTRL_U)) THEN
87 WHERE(ICNTRL_U(:) >= 0) ICNTRL(1:20) = ICNTRL_U(:)
88 ENDIF
89 IF (PRESENT(RCNTRL_U)) THEN
90 WHERE(RCNTRL_U(:) >= 0) RCNTRL(1:20) = RCNTRL_U(:)
91 ENDIF
92
93
94 CALL RosenbrockADJ(Y, NADJ, Lambda, &
95 TIN,TOUT, &
96 ATOL,RTOL, &
97 RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
98
99
100 ! N_stp = N_stp + ICNTRL(istp)
101 ! N_acc = N_acc + ICNTRL(iacc)
102 ! N_rej = N_rej + ICNTRL(irej)
103 ! N_sng = N_sng + ICNTRL(isng)
104 ! PRINT*,'Step=',N_stp,' Acc=',N_acc,' Rej=',N_rej, &
105 ! ' Singular=',N_sng
106
107 IF (IERR < 0) THEN
108 print *,'RosenbrockADJ: Unsucessful step at T=', &
109 TIN,' (IERR=',IERR,')'
110 ENDIF
111
112 STEPMIN = RCNTRL(ihexit)
113 ! if optional parameters are given for output
114 ! copy to them to return information
115 IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(1:20)
116 IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(1:20)
117
118 END SUBROUTINE INTEGRATE_ADJ
119
120 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
121 SUBROUTINE ros_AllocateDBuffers( S )
122 !~~~> Allocate buffer space for discrete adjoint
123 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
124 INTEGER :: i, S
125
126 ALLOCATE( buf_H(bufsize), STAT=i )
127 IF (i/=0) THEN
128 PRINT*,'Failed allocation of buffer H'; STOP
129 END IF
130 ALLOCATE( buf_T(bufsize), STAT=i )
131 IF (i/=0) THEN
132 PRINT*,'Failed allocation of buffer T'; STOP
133 END IF
134 ALLOCATE( buf_Y(NVAR*S,bufsize), STAT=i )
135 IF (i/=0) THEN
136 PRINT*,'Failed allocation of buffer Y'; STOP
137 END IF
138 ALLOCATE( buf_K(NVAR*S,bufsize), STAT=i )
139 IF (i/=0) THEN
140 PRINT*,'Failed allocation of buffer K'; STOP
141 END IF
142 ALLOCATE( buf_J(LU_NONZERO,bufsize), STAT=i )
143 IF (i/=0) THEN
144 PRINT*,'Failed allocation of buffer J'; STOP
145 END IF
146
147 END SUBROUTINE ros_AllocateDBuffers
148
149
150 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
151 SUBROUTINE ros_FreeDBuffers
152 !~~~> Dallocate buffer space for discrete adjoint
153 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
154 INTEGER :: i, S
155
156 DEALLOCATE( buf_H, STAT=i )
157 IF (i/=0) THEN
158 PRINT*,'Failed deallocation of buffer H'; STOP
159 END IF
160 DEALLOCATE( buf_T, STAT=i )
161 IF (i/=0) THEN
162 PRINT*,'Failed deallocation of buffer T'; STOP
163 END IF
164 DEALLOCATE( buf_Y, STAT=i )
165 IF (i/=0) THEN
166 PRINT*,'Failed deallocation of buffer Y'; STOP
167 END IF
168 DEALLOCATE( buf_K, STAT=i )
169 IF (i/=0) THEN
170 PRINT*,'Failed deallocation of buffer K'; STOP
171 END IF
172 DEALLOCATE( buf_J, STAT=i )
173 IF (i/=0) THEN
174 PRINT*,'Failed deallocation of buffer J'; STOP
175 END IF
176
177 END SUBROUTINE ros_FreeDBuffers
178
179
180 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
181 SUBROUTINE ros_AllocateCBuffers
182 !~~~> Allocate buffer space for continuous adjoint
183 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
184 INTEGER :: i, S
185
186 ALLOCATE( buf_H(bufsize), STAT=i )
187 IF (i/=0) THEN
188 PRINT*,'Failed allocation of buffer H'; STOP
189 END IF
190 ALLOCATE( buf_T(bufsize), STAT=i )
191 IF (i/=0) THEN
192 PRINT*,'Failed allocation of buffer T'; STOP
193 END IF
194 ALLOCATE( buf_Y(NVAR,bufsize), STAT=i )
195 IF (i/=0) THEN
196 PRINT*,'Failed allocation of buffer Y'; STOP
197 END IF
198 ALLOCATE( buf_dY(NVAR,bufsize), STAT=i )
199 IF (i/=0) THEN
200 PRINT*,'Failed allocation of buffer dY'; STOP
201 END IF
202 ALLOCATE( buf_d2Y(NVAR,bufsize), STAT=i )
203 IF (i/=0) THEN
204 PRINT*,'Failed allocation of buffer d2Y'; STOP
205 END IF
206
207 END SUBROUTINE ros_AllocateCBuffers
208
209
210 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
211 SUBROUTINE ros_FreeCBuffers
212 !~~~> Dallocate buffer space for continuous adjoint
213 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
214 INTEGER :: i, S
215
216 DEALLOCATE( buf_H, STAT=i )
217 IF (i/=0) THEN
218 PRINT*,'Failed deallocation of buffer H'; STOP
219 END IF
220 DEALLOCATE( buf_T, STAT=i )
221 IF (i/=0) THEN
222 PRINT*,'Failed deallocation of buffer T'; STOP
223 END IF
224 DEALLOCATE( buf_Y, STAT=i )
225 IF (i/=0) THEN
226 PRINT*,'Failed deallocation of buffer Y'; STOP
227 END IF
228 DEALLOCATE( buf_dY, STAT=i )
229 IF (i/=0) THEN
230 PRINT*,'Failed deallocation of buffer dY'; STOP
231 END IF
232 DEALLOCATE( buf_d2Y, STAT=i )
233 IF (i/=0) THEN
234 PRINT*,'Failed deallocation of buffer d2Y'; STOP
235 END IF
236
237 END SUBROUTINE ros_FreeCBuffers
238
239 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
240 SUBROUTINE ros_DPush( S, T, H, Ystage, K )!, Jcb )
241 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
242 !~~~> Saves the next trajectory snapshot for discrete adjoints
243 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
244 INTEGER :: S ! no of stages
245 KPP_REAL :: T, H, Ystage(NVAR*S), K(NVAR*S) !, Jcb(LU_NONZERO)
246
247 stack_ptr = stack_ptr + 1
248 IF ( stack_ptr > bufsize ) THEN
249 PRINT*,'Push failed: buffer overflow'
250 STOP
251 END IF
252 buf_H( stack_ptr ) = H
253 buf_T( stack_ptr ) = T
254 CALL WCOPY(NVAR*S,Ystage,1,buf_Y(1,stack_ptr),1)
255 CALL WCOPY(NVAR*S,K,1,buf_K(1,stack_ptr),1)
256 !CALL WCOPY(LU_NONZERO,Jcb,1,buf_J(1,stack_ptr),1)
257
258 END SUBROUTINE ros_DPush
259
260
261 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
262 SUBROUTINE ros_DPop( S, T, H, Ystage, K ) !, Jcb )
263 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
264 !~~~> Retrieves the next trajectory snapshot for discrete adjoints
265 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
266
267 INTEGER :: S ! no of stages
268 KPP_REAL :: T, H, Ystage(NVAR*S), K(NVAR*S) ! , Jcb(LU_NONZERO)
269
270 IF ( stack_ptr <= 0 ) THEN
271 PRINT*,'Pop failed: empty buffer'
272 STOP
273 END IF
274 H = buf_H( stack_ptr )
275 T = buf_T( stack_ptr )
276 CALL WCOPY(NVAR*S,buf_Y(1,stack_ptr),1,Ystage,1)
277 CALL WCOPY(NVAR*S,buf_K(1,stack_ptr),1,K,1)
278 !CALL WCOPY(LU_NONZERO,buf_J(1,stack_ptr),1,Jcb,1)
279
280 stack_ptr = stack_ptr - 1
281
282 END SUBROUTINE ros_DPop
283
284 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
285 SUBROUTINE ros_CPush( T, H, Y, dY, d2Y )
286 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
287 !~~~> Saves the next trajectory snapshot for discrete adjoints
288 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
289
290 INTEGER :: S ! no of stages
291 KPP_REAL :: T, H, Y(NVAR), dY(NVAR), d2Y(NVAR)
292
293 stack_ptr = stack_ptr + 1
294 IF ( stack_ptr > bufsize ) THEN
295 PRINT*,'Push failed: buffer overflow'
296 STOP
297 END IF
298 buf_H( stack_ptr ) = H
299 buf_T( stack_ptr ) = T
300 CALL WCOPY(NVAR,Y,1,buf_Y(1,stack_ptr),1)
301 CALL WCOPY(NVAR,dY,1,buf_dY(1,stack_ptr),1)
302 CALL WCOPY(NVAR,d2Y,1,buf_d2Y(1,stack_ptr),1)
303
304 END SUBROUTINE ros_CPush
305
306
307 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
308 SUBROUTINE ros_CPop( T, H, Y, dY, d2Y )
309 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
310 !~~~> Retrieves the next trajectory snapshot for discrete adjoints
311 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
312
313 INTEGER :: S ! no of stages
314 KPP_REAL :: T, H, Y(NVAR), dY(NVAR), d2Y(NVAR)
315
316 IF ( stack_ptr <= 0 ) THEN
317 PRINT*,'Pop failed: empty buffer'
318 STOP
319 END IF
320 H = buf_H( stack_ptr )
321 T = buf_T( stack_ptr )
322 CALL WCOPY(NVAR,buf_Y(1,stack_ptr),1,Y,1)
323 CALL WCOPY(NVAR,buf_dY(1,stack_ptr),1,dY,1)
324 CALL WCOPY(NVAR,buf_d2Y(1,stack_ptr),1,d2Y,1)
325
326 stack_ptr = stack_ptr - 1
327
328 END SUBROUTINE ros_CPop
329
330
331
332 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
333 SUBROUTINE RosenbrockADJ( Y, NADJ, Lambda, &
334 Tstart,Tend, &
335 AbsTol,RelTol, &
336 RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
337 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
338 !
339 ! ADJ = Adjoint of the Tangent Linear Model of a RosenbrockADJ Method
340 !
341 ! Solves the system y'=F(t,y) using a RosenbrockADJ method defined by:
342 !
343 ! G = 1/(H*gamma(1)) - Jac(t0,Y0)
344 ! T_i = t0 + Alpha(i)*H
345 ! Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
346 ! G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
347 ! gamma(i)*dF/dT(t0, Y0)
348 ! Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
349 !
350 ! For details on RosenbrockADJ methods and their implementation consult:
351 ! E. Hairer and G. Wanner
352 ! "Solving ODEs II. Stiff and differential-algebraic problems".
353 ! Springer series in computational mathematics, Springer-Verlag, 1996.
354 ! The codes contained in the book inspired this implementation.
355 !
356 ! (C) Adrian Sandu, August 2004
357 ! Virginia Polytechnic Institute and State University
358 ! Contact: sandu@cs.vt.edu
359 ! This implementation is part of KPP - the Kinetic PreProcessor
360 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
361 !
362 !~~~> INPUT ARGUMENTS:
363 !
364 !- Y(NVAR) = vector of initial conditions (at T=Tstart)
365 ! NADJ -> dimension of linearized system,
366 ! i.e. the number of sensitivity coefficients
367 !- Lambda(NVAR,NADJ) -> vector of initial sensitivity conditions (at T=Tstart)
368 !- [Tstart,Tend] = time range of integration
369 ! (if Tstart>Tend the integration is performed backwards in time)
370 !- RelTol, AbsTol = user precribed accuracy
371 !- SUBROUTINE Fun( T, Y, Ydot ) = ODE function,
372 ! returns Ydot = Y' = F(T,Y)
373 !- SUBROUTINE Jac( T, Y, Jcb ) = Jacobian of the ODE function,
374 ! returns Jcb = dF/dY
375 !- ICNTRL(1:10) = integer inputs parameters
376 !- RCNTRL(1:10) = real inputs parameters
377 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
378 !
379 !~~~> OUTPUT ARGUMENTS:
380 !
381 !- Y(NVAR) -> vector of final states (at T->Tend)
382 !- Lambda(NVAR,NADJ) -> vector of final sensitivities (at T=Tend)
383 !- ICNTRL(11:20) -> integer output parameters
384 !- RCNTRL(11:20) -> real output parameters
385 !- IERR -> job status upon return
386 ! - succes (positive value) or failure (negative value) -
387 ! = 1 : Success
388 ! = -1 : Improper value for maximal no of steps
389 ! = -2 : Selected RosenbrockADJ method not implemented
390 ! = -3 : Hmin/Hmax/Hstart must be positive
391 ! = -4 : FacMin/FacMax/FacRej must be positive
392 ! = -5 : Improper tolerance values
393 ! = -6 : No of steps exceeds maximum bound
394 ! = -7 : Step size too small
395 ! = -8 : Matrix is repeatedly singular
396 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
397 !
398 !~~~> INPUT PARAMETERS:
399 !
400 ! Note: For input parameters equal to zero the default values of the
401 ! corresponding variables are used.
402 !
403 ! ICNTRL(1) = 1: F = F(y) Independent of T (AUTONOMOUS)
404 ! = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
405 ! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
406 ! = 1: AbsTol, RelTol are scalars
407 ! ICNTRL(3) -> maximum number of integration steps
408 ! For ICNTRL(3)=0) the default value of 100000 is used
409 !
410 ! ICNTRL(4) -> selection of a particular Rosenbrock method
411 ! = 0 : default method is Rodas3
412 ! = 1 : method is Ros2
413 ! = 2 : method is Ros3
414 ! = 3 : method is Ros4
415 ! = 4 : method is Rodas3
416 ! = 5: method is Rodas4
417 !
418 ! ICNTRL(5) -> Type of adjoint algorithm
419 ! = 0 : default is discrete adjoint ( of method ICNTRL(4) )
420 ! = 1 : no adjoint
421 ! = 2 : discrete adjoint ( of method ICNTRL(4) )
422 ! = 3 : fully adaptive continuous adjoint ( with method ICNTRL(6) )
423 ! = 4 : simplified continuous adjoint ( with method ICNTRL(6) )
424 !
425 ! ICNTRL(6) -> selection of a particular Rosenbrock method for the
426 ! continuous adjoint integration - for cts adjoint it
427 ! can be different than the forward method ICNTRL(4)
428 ! Note 1: to avoid interpolation errors (which can be huge!)
429 ! it is recommended to use only ICNTRL(6) = 1 or 4
430 ! Note 2: the performance of the full continuous adjoint
431 ! strongly depends on the forward solution accuracy Abs/RelTol
432 !
433 ! RCNTRL(1) -> Hmin, lower bound for the integration step size
434 ! It is strongly recommended to keep Hmin = ZERO
435 ! RCNTRL(2) -> Hmax, upper bound for the integration step size
436 ! RCNTRL(3) -> Hstart, starting value for the integration step size
437 !
438 ! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2)
439 ! RCNTRL(5) -> FacMin,upper bound on step increase factor (default=6)
440 ! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections
441 ! (default=0.1)
442 ! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller
443 ! than the predicted value (default=0.9)
444 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
445 !
446 !~~~> OUTPUT PARAMETERS:
447 !
448 ! Note: each call to RosenbrockADJ adds the corrent no. of fcn calls
449 ! to previous value of ISTATUS(1), and similar for the other params.
450 ! Set ISTATUS(1:10) = 0 before call to avoid this accumulation.
451 !
452 ! ISTATUS(1) = No. of function calls
453 ! ISTATUS(2) = No. of jacobian calls
454 ! ISTATUS(3) = No. of steps
455 ! ISTATUS(4) = No. of accepted steps
456 ! ISTATUS(5) = No. of rejected steps (except at the beginning)
457 ! ISTATUS(6) = No. of LU decompositions
458 ! ISTATUS(7) = No. of forward/backward substitutions
459 ! ISTATUS(8) = No. of singular matrix decompositions
460 !
461 ! RSTATUS(1) -> Texit, the time corresponding to the
462 ! computed Y upon return
463 ! RSTATUS(2) -> Hexit, last accepted step before exit
464 ! For multiple restarts, use Hexit as Hstart in the following run
465 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
466
467 IMPLICIT NONE
468
469 !~~~> Arguments
470 KPP_REAL, INTENT(INOUT) :: Y(NVAR)
471 INTEGER, INTENT(IN) :: NADJ
472 KPP_REAL, INTENT(INOUT) :: Lambda(NVAR,NADJ)
473 KPP_REAL, INTENT(IN) :: Tstart,Tend
474 KPP_REAL, INTENT(IN) :: AbsTol(NVAR),RelTol(NVAR)
475 INTEGER, INTENT(IN) :: ICNTRL(10)
476 KPP_REAL, INTENT(IN) :: RCNTRL(10)
477 INTEGER, INTENT(INOUT) :: ISTATUS(10)
478 KPP_REAL, INTENT(INOUT) :: RSTATUS(10)
479 INTEGER, INTENT(OUT) :: IERR
480 !~~~> The method parameters
481 INTEGER, PARAMETER :: Smax = 6
482 INTEGER :: Method, ros_S
483 KPP_REAL, DIMENSION(Smax) :: ros_M, ros_E, ros_Alpha, ros_Gamma
484 KPP_REAL, DIMENSION(Smax*(Smax-1)/2) :: ros_A, ros_C
485 KPP_REAL :: ros_ELO
486 LOGICAL, DIMENSION(Smax) :: ros_NewF
487 CHARACTER(LEN=12) :: ros_Name
488 !~~~> Local variables
489 KPP_REAL :: Roundoff, FacMin, FacMax, FacRej, FacSafe
490 KPP_REAL :: Hmin, Hmax, Hstart, Hexit
491 KPP_REAL :: Texit
492 INTEGER :: i, UplimTol, Max_no_steps
493 INTEGER :: AdjointType, CadjMethod
494 LOGICAL :: Autonomous, VectorTol
495 !~~~> Parameters
496 KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
497 KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
498
499 !~~~> Initialize statistics
500 Nfun = ISTATUS(ifun)
501 Njac = ISTATUS(ijac)
502 Nstp = ISTATUS(istp)
503 Nacc = ISTATUS(iacc)
504 Nrej = ISTATUS(irej)
505 Ndec = ISTATUS(idec)
506 Nsol = ISTATUS(isol)
507 Nsng = ISTATUS(isng)
508
509 !~~~> Autonomous or time dependent ODE. Default is time dependent.
510 Autonomous = .NOT.(ICNTRL(1) == 0)
511
512 !~~~> For Scalar tolerances (ICNTRL(2).NE.0) the code uses AbsTol(1) and RelTol(1)
513 ! For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
514 IF (ICNTRL(2) == 0) THEN
515 VectorTol = .TRUE.
516 UplimTol = NVAR
517 ELSE
518 VectorTol = .FALSE.
519 UplimTol = 1
520 END IF
521
522 !~~~> The maximum number of steps admitted
523 IF (ICNTRL(3) == 0) THEN
524 Max_no_steps = bufsize - 1
525 ELSEIF (Max_no_steps > 0) THEN
526 Max_no_steps=ICNTRL(3)
527 ELSE
528 PRINT * ,'User-selected max no. of steps: ICNTRL(3)=',ICNTRL(3)
529 CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR)
530 RETURN
531 END IF
532
533 !~~~> The particular Rosenbrock method chosen
534 IF (ICNTRL(4) == 0) THEN
535 Method = 5
536 ELSEIF ( (ICNTRL(4) >= 1).AND.(ICNTRL(4) <= 5) ) THEN
537 Method = ICNTRL(4)
538 ELSE
539 PRINT * , 'User-selected Rosenbrock method: ICNTRL(4)=', Method
540 CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
541 RETURN
542 END IF
543
544 !~~~> Discrete or continuous adjoint formulation
545 IF ( ICNTRL(5) == 0 ) THEN
546 AdjointType = Adj_discrete
547 ELSEIF ( (ICNTRL(5) >= 1).AND.(ICNTRL(5) <= 4) ) THEN
548 AdjointType = ICNTRL(5)
549 ELSE
550 PRINT * , 'User-selected adjoint type: ICNTRL(5)=', AdjointType
551 CALL ros_ErrorMsg(-9,Tstart,ZERO,IERR)
552 RETURN
553 END IF
554
555 !~~~> The particular Rosenbrock method chosen for integrating the cts adjoint
556 IF (ICNTRL(6) == 0) THEN
557 CadjMethod = 4
558 ELSEIF ( (ICNTRL(6) >= 1).AND.(ICNTRL(6) <= 5) ) THEN
559 CadjMethod = ICNTRL(4)
560 ELSE
561 PRINT * , 'User-selected CADJ Rosenbrock method: ICNTRL(6)=', Method
562 CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
563 RETURN
564 END IF
565
566
567 !~~~> Unit roundoff (1+Roundoff>1)
568 Roundoff = WLAMCH('E')
569
570 !~~~> Lower bound on the step size: (positive value)
571 IF (RCNTRL(1) == ZERO) THEN
572 Hmin = ZERO
573 ELSEIF (RCNTRL(1) > ZERO) THEN
574 Hmin = RCNTRL(1)
575 ELSE
576 PRINT * , 'User-selected Hmin: RCNTRL(1)=', RCNTRL(1)
577 CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
578 RETURN
579 END IF
580 !~~~> Upper bound on the step size: (positive value)
581 IF (RCNTRL(2) == ZERO) THEN
582 Hmax = ABS(Tend-Tstart)
583 ELSEIF (RCNTRL(2) > ZERO) THEN
584 Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart))
585 ELSE
586 PRINT * , 'User-selected Hmax: RCNTRL(2)=', RCNTRL(2)
587 CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
588 RETURN
589 END IF
590 !~~~> Starting step size: (positive value)
591 IF (RCNTRL(3) == ZERO) THEN
592 Hstart = MAX(Hmin,DeltaMin)
593 ELSEIF (RCNTRL(3) > ZERO) THEN
594 Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart))
595 ELSE
596 PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
597 CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
598 RETURN
599 END IF
600 !~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax
601 IF (RCNTRL(4) == ZERO) THEN
602 FacMin = 0.2d0
603 ELSEIF (RCNTRL(4) > ZERO) THEN
604 FacMin = RCNTRL(4)
605 ELSE
606 PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
607 CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
608 RETURN
609 END IF
610 IF (RCNTRL(5) == ZERO) THEN
611 FacMax = 6.0d0
612 ELSEIF (RCNTRL(5) > ZERO) THEN
613 FacMax = RCNTRL(5)
614 ELSE
615 PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
616 CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
617 RETURN
618 END IF
619 !~~~> FacRej: Factor to decrease step after 2 succesive rejections
620 IF (RCNTRL(6) == ZERO) THEN
621 FacRej = 0.1d0
622 ELSEIF (RCNTRL(6) > ZERO) THEN
623 FacRej = RCNTRL(6)
624 ELSE
625 PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
626 CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
627 RETURN
628 END IF
629 !~~~> FacSafe: Safety Factor in the computation of new step size
630 IF (RCNTRL(7) == ZERO) THEN
631 FacSafe = 0.9d0
632 ELSEIF (RCNTRL(7) > ZERO) THEN
633 FacSafe = RCNTRL(7)
634 ELSE
635 PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
636 CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
637 RETURN
638 END IF
639 !~~~> Check if tolerances are reasonable
640 DO i=1,UplimTol
641 IF ( (AbsTol(i) <= ZERO) .OR. (RelTol(i) <= 10.d0*Roundoff) &
642 .OR. (RelTol(i) >= 1.0d0) ) THEN
643 PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
644 PRINT * , ' RelTol(',i,') = ',RelTol(i)
645 CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR)
646 RETURN
647 END IF
648 END DO
649
650
651 !~~~> Initialize the particular RosenbrockADJ method
652 SELECT CASE (Method)
653 CASE (1)
654 CALL Ros2(ros_S, ros_A, ros_C, ros_M, ros_E, &
655 ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
656 CASE (2)
657 CALL Ros3(ros_S, ros_A, ros_C, ros_M, ros_E, &
658 ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
659 CASE (3)
660 CALL Ros4(ros_S, ros_A, ros_C, ros_M, ros_E, &
661 ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
662 CASE (4)
663 CALL Rodas3(ros_S, ros_A, ros_C, ros_M, ros_E, &
664 ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
665 CASE (5)
666 CALL Rodas4(ros_S, ros_A, ros_C, ros_M, ros_E, &
667 ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
668 CASE DEFAULT
669 PRINT * , 'Unknown Rosenbrock method: ICNTRL(4)=', Method
670 CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
671 RETURN
672 END SELECT
673
674 !~~~> Allocate checkpoint space or open checkpoint files
675 IF (AdjointType == Adj_discrete) THEN
676 CALL ros_AllocateDBuffers( ros_S )
677 ELSEIF ( (AdjointType == Adj_continuous).OR. &
678 (AdjointType == Adj_simple_continuous) ) THEN
679 CALL ros_AllocateCBuffers
680 END IF
681
682 !~~~> CALL Forward Rosenbrock method
683 CALL ros_FwdInt(Y,Tstart,Tend,Texit, &
684 AbsTol, RelTol, &
685 ! RosenbrockADJ method coefficients
686 ros_S, ros_M, ros_E, ros_A, ros_C, &
687 ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
688 ! Integration parameters
689 Autonomous, VectorTol, AdjointType, &
690 Max_no_steps, &
691 Roundoff, Hmin, Hmax, Hstart, Hexit, &
692 FacMin, FacMax, FacRej, FacSafe, &
693 ! Error indicator
694 IERR)
695
696 PRINT*,'FORWARD STATISTICS'
697 PRINT*,'Step=',Nstp,' Acc=',Nacc, &
698 ' Rej=',Nrej, ' Singular=',Nsng
699 Nstp = 0
700 Nacc = 0
701 Nrej = 0
702 Nsng = 0
703
704 !~~~> If Forward integration failed return
705 IF (IERR<0) RETURN
706
707 !~~~> Initialize the particular Rosenbrock method for continuous adjoint
708 IF ( (AdjointType == Adj_continuous).OR. &
709 (AdjointType == Adj_simple_continuous) ) THEN
710 SELECT CASE (CadjMethod)
711 CASE (1)
712 CALL Ros2(ros_S, ros_A, ros_C, ros_M, ros_E, &
713 ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
714 CASE (2)
715 CALL Ros3(ros_S, ros_A, ros_C, ros_M, ros_E, &
716 ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
717 CASE (3)
718 CALL Ros4(ros_S, ros_A, ros_C, ros_M, ros_E, &
719 ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
720 CASE (4)
721 CALL Rodas3(ros_S, ros_A, ros_C, ros_M, ros_E, &
722 ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
723 CASE (5)
724 CALL Rodas4(ros_S, ros_A, ros_C, ros_M, ros_E, &
725 ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
726 CASE DEFAULT
727 PRINT * , 'Unknown Rosenbrock method: ICNTRL(4)=', Method
728 CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
729 RETURN
730 END SELECT
731 END IF
732
733 SELECT CASE (AdjointType)
734 CASE (Adj_discrete)
735 CALL ros_DadjInt ( &
736 NADJ, Lambda, &
737 Tstart, Tend, Texit, &
738 AbsTol, RelTol, &
739 ros_S, ros_M, ros_E, ros_A, ros_C, &
740 ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
741 Autonomous, VectorTol, Max_no_steps, &
742 Roundoff, Hmin, Hmax, Hstart, &
743 FacMin, FacMax, FacRej, FacSafe, &
744 IERR )
745 CASE (Adj_continuous)
746 CALL ros_CadjInt ( &
747 NADJ, Lambda, &
748 Tend, Tstart, Texit, &
749 AbsTol, RelTol, &
750 ros_S, ros_M, ros_E, ros_A, ros_C, &
751 ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
752 Autonomous, VectorTol, AdjointType, &
753 100000, &
754 Roundoff, Hmin, Hmax, Hstart, Hexit, &
755 FacMin, FacMax, FacRej, FacSafe, &
756 IERR )
757 CASE (Adj_simple_continuous)
758 CALL ros_SimpleCadjInt ( &
759 NADJ, Lambda, &
760 Tstart, Tend, Texit, &
761 AbsTol, RelTol, &
762 ros_S, ros_M, ros_E, ros_A, ros_C, &
763 ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
764 Autonomous, VectorTol, AdjointType, &
765 Max_no_steps, &
766 Roundoff, Hmin, Hmax, Hstart, &
767 FacMin, FacMax, FacRej, FacSafe, &
768 IERR )
769 END SELECT ! AdjointType
770
771 PRINT*,'ADJOINT STATISTICS'
772 PRINT*,'Step=',Nstp,' Acc=',Nacc, &
773 ' Rej=',Nrej, ' Singular=',Nsng
774
775 !~~~> Free checkpoint space or close checkpoint files
776 IF (AdjointType == Adj_discrete) THEN
777 CALL ros_FreeDBuffers
778 ELSEIF ( (AdjointType == Adj_continuous) .OR. &
779 (AdjointType == Adj_simple_continuous) ) THEN
780 CALL ros_FreeCBuffers
781 END IF
782
783 !~~~> Collect run statistics
784 ISTATUS(ifun) = Nfun
785 ISTATUS(ijac) = Njac
786 ISTATUS(istp) = Nstp
787 ISTATUS(iacc) = Nacc
788 ISTATUS(irej) = Nrej
789 ISTATUS(idec) = Ndec
790 ISTATUS(isol) = Nsol
791 ISTATUS(isng) = Nsng
792 !~~~> Last T and H
793 RSTATUS(itexit) = Texit
794 RSTATUS(ihexit) = Hexit
795
796
797 END SUBROUTINE RosenbrockADJ
798 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
799
800
801 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
802 SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
803 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
804 ! Handles all error messages
805 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
806
807 KPP_REAL, INTENT(IN) :: T, H
808 INTEGER, INTENT(IN) :: Code
809 INTEGER, INTENT(OUT) :: IERR
810
811 IERR = Code
812 PRINT * , &
813 'Forced exit from RosenbrockADJ due to the following error:'
814
815 SELECT CASE (Code)
816 CASE (-1)
817 PRINT * , '--> Improper value for maximal no of steps'
818 CASE (-2)
819 PRINT * , '--> Selected RosenbrockADJ method not implemented'
820 CASE (-3)
821 PRINT * , '--> Hmin/Hmax/Hstart must be positive'
822 CASE (-4)
823 PRINT * , '--> FacMin/FacMax/FacRej must be positive'
824 CASE (-5)
825 PRINT * , '--> Improper tolerance values'
826 CASE (-6)
827 PRINT * , '--> No of steps exceeds maximum buffer bound'
828 CASE (-7)
829 PRINT * , '--> Step size too small: T + 10*H = T', &
830 ' or H < Roundoff'
831 CASE (-8)
832 PRINT * , '--> Matrix is repeatedly singular'
833 CASE (-9)
834 PRINT * , '--> Improper type of adjoint selected'
835 CASE DEFAULT
836 PRINT *, 'Unknown Error code: ', Code
837 END SELECT
838
839 PRINT *, "T=", T, "and H=", H
840
841 END SUBROUTINE ros_ErrorMsg
842
843
844
845 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
846 SUBROUTINE ros_FwdInt (Y, &
847 Tstart, Tend, T, &
848 AbsTol, RelTol, &
849 !~~~> RosenbrockADJ method coefficients
850 ros_S, ros_M, ros_E, ros_A, ros_C, &
851 ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
852 !~~~> Integration parameters
853 Autonomous, VectorTol, AdjointType, &
854 Max_no_steps, &
855 Roundoff, Hmin, Hmax, Hstart, Hexit, &
856 FacMin, FacMax, FacRej, FacSafe, &
857 !~~~> Error indicator
858 IERR )
859 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
860 ! Template for the implementation of a generic RosenbrockADJ method
861 ! defined by ros_S (no of stages)
862 ! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
863 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
864
865 IMPLICIT NONE
866
867 !~~~> Input: the initial condition at Tstart; Output: the solution at T
868 KPP_REAL, INTENT(INOUT) :: Y(NVAR)
869 !~~~> Input: integration interval
870 KPP_REAL, INTENT(IN) :: Tstart,Tend
871 !~~~> Output: time at which the solution is returned (T=Tend if success)
872 KPP_REAL, INTENT(OUT) :: T
873 !~~~> Input: tolerances
874 KPP_REAL, INTENT(IN) :: AbsTol(NVAR), RelTol(NVAR)
875 !~~~> Input: The RosenbrockADJ method parameters
876 INTEGER, INTENT(IN) :: ros_S
877 KPP_REAL, INTENT(IN) :: ros_M(ros_S), ros_E(ros_S), &
878 ros_Alpha(ros_S), ros_A(ros_S*(ros_S-1)/2), &
879 ros_Gamma(ros_S), ros_C(ros_S*(ros_S-1)/2), ros_ELO
880 LOGICAL, INTENT(IN) :: ros_NewF(ros_S)
881 !~~~> Input: integration parameters
882 LOGICAL, INTENT(IN) :: Autonomous, VectorTol
883 INTEGER, INTENT(IN) :: AdjointType
884 KPP_REAL, INTENT(IN) :: Hstart, Hmin, Hmax
885 INTEGER, INTENT(IN) :: Max_no_steps
886 KPP_REAL, INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
887 !~~~> Output: last accepted step
888 KPP_REAL, INTENT(OUT) :: Hexit
889 !~~~> Output: Error indicator
890 INTEGER, INTENT(OUT) :: IERR
891 ! ~~~~ Local variables
892 KPP_REAL :: Ynew(NVAR), Fcn0(NVAR), Fcn(NVAR)
893 KPP_REAL :: K(NVAR*ros_S), dFdT(NVAR)
894 KPP_REAL, DIMENSION(:), POINTER :: Ystage
895 KPP_REAL :: Jac0(LU_NONZERO), Ghimj(LU_NONZERO)
896 KPP_REAL :: H, Hnew, HC, HG, Fac, Tau
897 KPP_REAL :: Err, Yerr(NVAR)
898 INTEGER :: Pivot(NVAR), Direction, ioffset, i, j, istage
899 LOGICAL :: RejectLastH, RejectMoreH, Singular
900 !~~~> Local parameters
901 KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
902 !~~~> Locally called functions
903 ! KPP_REAL WLAMCH
904 ! EXTERNAL WLAMCH
905 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
906
907 !~~~> Allocate stage vector buffer if needed
908 IF (AdjointType == Adj_discrete) THEN ! Save stage solution
909 ALLOCATE(Ystage(NVAR*ros_S), STAT=i)
910 IF (i/=0) THEN
911 PRINT*,'Allocation of Ystage failed'
912 STOP
913 END IF
914 END IF
915
916 !~~~> Initial preparations
917 T = Tstart
918 Hexit = 0.0_dp
919 H = MIN(Hstart,Hmax)
920 IF (ABS(H) <= 10.0_dp*Roundoff) H = DeltaMin
921
922 IF (Tend >= Tstart) THEN
923 Direction = +1
924 ELSE
925 Direction = -1
926 END IF
927
928 RejectLastH=.FALSE.
929 RejectMoreH=.FALSE.
930
931 !~~~> Time loop begins below
932
933 TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff <= ZERO) &
934 .OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) )
935
936 IF ( Nstp > Max_no_steps ) THEN ! Too many steps
937 CALL ros_ErrorMsg(-6,T,H,IERR)
938 RETURN
939 END IF
940 IF ( ((T+0.1d0*H) == T).OR.(H <= Roundoff) ) THEN ! Step size too small
941 CALL ros_ErrorMsg(-7,T,H,IERR)
942 RETURN
943 END IF
944
945 !~~~> Limit H if necessary to avoid going beyond Tend
946 Hexit = H
947 H = MIN(H,ABS(Tend-T))
948
949 !~~~> Compute the function at current time
950 CALL FunTemplate(T,Y,Fcn0)
951
952 !~~~> Compute the function derivative with respect to T
953 IF (.NOT.Autonomous) THEN
954 CALL ros_FunTimeDerivative ( T, Roundoff, Y, &
955 Fcn0, dFdT )
956 END IF
957
958 !~~~> Compute the Jacobian at current time
959 CALL JacTemplate(T,Y,Jac0)
960
961 !~~~> Repeat step calculation until current step accepted
962 UntilAccepted: DO
963
964 CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1), &
965 Jac0,Ghimj,Pivot,Singular)
966 IF (Singular) THEN ! More than 5 consecutive failed decompositions
967 CALL ros_ErrorMsg(-8,T,H,IERR)
968 RETURN
969 END IF
970
971 !~~~> Compute the stages
972 Stage: DO istage = 1, ros_S
973
974 ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR)
975 ioffset = NVAR*(istage-1)
976
977 ! For the 1st istage the function has been computed previously
978 IF ( istage == 1 ) THEN
979 CALL WCOPY(NVAR,Fcn0,1,Fcn,1)
980 IF (AdjointType == Adj_discrete) THEN ! Save stage solution
981 CALL WCOPY(NVAR,Y,1,Ystage(1),1)
982 END IF
983 ! istage>1 and a new function evaluation is needed at the current istage
984 ELSEIF ( ros_NewF(istage) ) THEN
985 CALL WCOPY(NVAR,Y,1,Ynew,1)
986 DO j = 1, istage-1
987 CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), &
988 K(NVAR*(j-1)+1),1,Ynew,1)
989 END DO
990 Tau = T + ros_Alpha(istage)*Direction*H
991 CALL FunTemplate(Tau,Ynew,Fcn)
992 IF (AdjointType == Adj_discrete) THEN ! Save stage solution
993 CALL WCOPY(NVAR,Ynew,1,Ystage(ioffset+1),1)
994 END IF
995 END IF ! if istage == 1 elseif ros_NewF(istage)
996 CALL WCOPY(NVAR,Fcn,1,K(ioffset+1),1)
997 DO j = 1, istage-1
998 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H)
999 CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1),1,K(ioffset+1),1)
1000 END DO
1001 IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
1002 HG = Direction*H*ros_Gamma(istage)
1003 CALL WAXPY(NVAR,HG,dFdT,1,K(ioffset+1),1)
1004 END IF
1005 CALL ros_Solve('N', Ghimj, Pivot, K(ioffset+1))
1006
1007 END DO Stage
1008
1009
1010 !~~~> Compute the new solution
1011 CALL WCOPY(NVAR,Y,1,Ynew,1)
1012 DO j=1,ros_S
1013 CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1),1,Ynew,1)
1014 END DO
1015
1016 !~~~> Compute the error estimation
1017 CALL WSCAL(NVAR,ZERO,Yerr,1)
1018 DO j=1,ros_S
1019 CALL WAXPY(NVAR,ros_E(j),K(NVAR*(j-1)+1),1,Yerr,1)
1020 END DO
1021 Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
1022
1023 !~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
1024 Fac = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
1025 Hnew = H*Fac
1026
1027 !~~~> Check the error magnitude and adjust step size
1028 Nstp = Nstp+1
1029 IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step
1030 Nacc = Nacc+1
1031 IF (AdjointType == Adj_discrete) THEN ! Save current state
1032 CALL ros_DPush( ros_S, T, H, Ystage, K ) !, Ghimj )
1033 ELSEIF ( (AdjointType == Adj_continuous) .OR. &
1034 (AdjointType == Adj_simple_continuous) ) THEN
1035 CALL Jac_SP_Vec( Jac0, Fcn0, K(1) )
1036 IF (.NOT. Autonomous) THEN
1037 CALL WAXPY(NVAR,ONE,dFdT,1,K(1),1)
1038 END IF
1039 CALL ros_CPush( T, H, Y, Fcn0, K(1) )
1040 END IF
1041 CALL WCOPY(NVAR,Ynew,1,Y,1)
1042 T = T + Direction*H
1043 Hnew = MAX(Hmin,MIN(Hnew,Hmax))
1044 IF (RejectLastH) THEN ! No step size increase after a rejected step
1045 Hnew = MIN(Hnew,H)
1046 END IF
1047 RejectLastH = .FALSE.
1048 RejectMoreH = .FALSE.
1049 H = Hnew
1050 EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
1051 ELSE !~~~> Reject step
1052 IF (RejectMoreH) THEN
1053 Hnew = H*FacRej
1054 END IF
1055 RejectMoreH = RejectLastH
1056 RejectLastH = .TRUE.
1057 H = Hnew
1058 IF (Nacc >= 1) THEN
1059 Nrej = Nrej+1
1060 END IF
1061 END IF ! Err <= 1
1062
1063 END DO UntilAccepted
1064
1065 END DO TimeLoop
1066
1067 !~~~> Save last state: only needed for continuous adjoint
1068 IF ( (AdjointType == Adj_continuous) .OR. &
1069 (AdjointType == Adj_simple_continuous) ) THEN
1070 CALL FunTemplate(T,Y,Fcn0)
1071 CALL JacTemplate(T,Y,Jac0)
1072 CALL Jac_SP_Vec( Jac0, Fcn0, K(1) )
1073 IF (.NOT. Autonomous) THEN
1074 CALL ros_FunTimeDerivative ( T, Roundoff, Y, &
1075 Fcn0, dFdT )
1076 CALL WAXPY(NVAR,ONE,dFdT,1,K(1),1)
1077 END IF
1078 CALL ros_CPush( T, H, Y, Fcn0, K(1) )
1079 !~~~> Deallocate stage buffer: only needed for discrete adjoint
1080 ELSEIF (AdjointType == Adj_discrete) THEN
1081 DEALLOCATE(Ystage, STAT=i)
1082 IF (i/=0) THEN
1083 PRINT*,'Deallocation of Ystage failed'
1084 STOP
1085 END IF
1086 END IF
1087
1088 !~~~> Succesful exit
1089 IERR = 1 !~~~> The integration was successful
1090
1091 PRINT*,'Nacc after fwd =',Nacc
1092
1093 END SUBROUTINE ros_FwdInt
1094
1095
1096
1097
1098
1099 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1100 SUBROUTINE ros_DadjInt ( &
1101 NADJ, Lambda, &
1102 Tstart, Tend, T, &
1103 AbsTol, RelTol, &
1104 !~~~> RosenbrockSOA method coefficients
1105 ros_S, ros_M, ros_E, ros_A, ros_C, &
1106 ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
1107 !~~~> Integration parameters
1108 Autonomous, VectorTol, Max_no_steps, &
1109 Roundoff, Hmin, Hmax, Hstart, &
1110 FacMin, FacMax, FacRej, FacSafe, &
1111 !~~~> Error indicator
1112 IERR )
1113 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1114 ! Template for the implementation of a generic RosenbrockSOA method
1115 ! defined by ros_S (no of stages)
1116 ! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1117 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1118
1119 IMPLICIT NONE
1120
1121 !~~~> Input: the initial condition at Tstart; Output: the solution at T
1122 INTEGER, INTENT(IN) :: NADJ
1123 !~~~> First order adjoint
1124 KPP_REAL, INTENT(INOUT) :: Lambda(NVAR,NADJ)
1125 !!~~~> Input: integration interval
1126 KPP_REAL, INTENT(IN) :: Tstart,Tend
1127 !~~~> Output: time at which the solution is returned (T=Tend if success)
1128 KPP_REAL, INTENT(OUT) :: T
1129 !~~~> Input: tolerances
1130 KPP_REAL, INTENT(IN) :: AbsTol(NVAR), RelTol(NVAR)
1131 !~~~> Input: The RosenbrockSOA method parameters
1132 INTEGER, INTENT(IN) :: ros_S
1133 KPP_REAL, INTENT(IN) :: ros_M(ros_S), ros_E(ros_S), &
1134 ros_Alpha(ros_S), ros_A(ros_S*(ros_S-1)/2), &
1135 ros_Gamma(ros_S), ros_C(ros_S*(ros_S-1)/2), ros_ELO
1136 LOGICAL, INTENT(IN) :: ros_NewF(ros_S)
1137 !~~~> Input: integration parameters
1138 LOGICAL, INTENT(IN) :: Autonomous, VectorTol
1139 KPP_REAL, INTENT(IN) :: Hstart, Hmin, Hmax
1140 INTEGER, INTENT(IN) :: Max_no_steps
1141 KPP_REAL, INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
1142 !~~~> Output: Error indicator
1143 INTEGER, INTENT(OUT) :: IERR
1144 ! ~~~~ Local variables
1145 KPP_REAL :: Ystage_adj(NVAR,NADJ)
1146 KPP_REAL :: dFdT(NVAR)
1147 KPP_REAL :: Ystage(NVAR*ros_S), K(NVAR*ros_S)
1148 KPP_REAL :: U(NVAR*ros_S,NADJ), V(NVAR*ros_S,NADJ)
1149 KPP_REAL :: Jac(LU_NONZERO), dJdT(LU_NONZERO), Ghimj(LU_NONZERO)
1150 KPP_REAL :: Hes0(NHESS), Hes1(NHESS), dHdT(NHESS)
1151 KPP_REAL :: Tmp(NVAR), Tmp2(NVAR)
1152 KPP_REAL :: H, HC, HA, Tau
1153 INTEGER :: Pivot(NVAR), Direction
1154 INTEGER :: i, j, m, istage, istart, jstart
1155 !~~~> Local parameters
1156 KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
1157 KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
1158 !~~~> Locally called functions
1159 ! KPP_REAL WLAMCH
1160 ! EXTERNAL WLAMCH
1161 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1162
1163
1164
1165 IF (Tend >= Tstart) THEN
1166 Direction = +1
1167 ELSE
1168 Direction = -1
1169 END IF
1170
1171 OPEN(55,file='KPP_ROOT_dadj.dat')
1172
1173 !~~~> Time loop begins below
1174 TimeLoop: DO WHILE ( stack_ptr > 0 )
1175
1176 !~~~> Recover checkpoints for stage values and vectors
1177 CALL ros_DPop( ros_S, T, H, Ystage, K ) !, Ghimj )
1178
1179 Nstp = Nstp+1
1180
1181 !~~~> Compute LU decomposition
1182 CALL JacTemplate(T,Ystage(1),Ghimj)
1183 CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
1184 Tau = ONE/(Direction*H*ros_Gamma(1))
1185 DO i=1,NVAR
1186 Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+Tau
1187 END DO
1188 CALL ros_Decomp( Ghimj, Pivot, j )
1189
1190 !~~~> Compute Hessian at the beginning of the interval
1191 CALL HessTemplate(T,Ystage(1),Hes0)
1192
1193 !~~~> Compute the stages
1194 Stage: DO istage = ros_S, 1, -1
1195
1196 !~~~> Current istage first entry
1197 istart = NVAR*(istage-1) + 1
1198
1199 !~~~> Compute U
1200 DO m = 1,NADJ
1201 CALL WCOPY(NVAR,Lambda(1,m),1,U(istart,m),1)
1202 CALL WSCAL(NVAR,ros_M(istage),U(istart,m),1)
1203 END DO ! m=1:NADJ
1204 DO j = istage+1, ros_S
1205 jstart = NVAR*(j-1) + 1
1206 HA = ros_A((j-1)*(j-2)/2+istage)
1207 HC = ros_C((j-1)*(j-2)/2+istage)/(Direction*H)
1208 DO m = 1,NADJ
1209 CALL WAXPY(NVAR,HA,V(jstart,m),1,U(istart,m),1)
1210 CALL WAXPY(NVAR,HC,U(jstart,m),1,U(istart,m),1)
1211 END DO ! m=1:NADJ
1212 END DO
1213 DO m = 1,NADJ
1214 CALL ros_Solve('T', Ghimj, Pivot, U(istart,m))
1215 END DO ! m=1:NADJ
1216 !~~~> Compute V
1217 Tau = T + ros_Alpha(istage)*Direction*H
1218 CALL JacTemplate(Tau,Ystage(istart),Jac)
1219 DO m = 1,NADJ
1220 CALL JacTR_SP_Vec(Jac,U(istart,m),V(istart,m))
1221 END DO ! m=1:NADJ
1222
1223 END DO Stage
1224
1225 IF (.NOT.Autonomous) THEN
1226 !~~~> Compute the Jacobian derivative with respect to T.
1227 ! Last "Jac" computed for stage 1
1228 CALL ros_JacTimeDerivative ( T, Roundoff, Ystage(1), &
1229 Jac, dJdT )
1230 END IF
1231
1232 !~~~> Compute the new solution
1233
1234 !~~~> Compute Lambda
1235 DO istage=1,ros_S
1236 istart = NVAR*(istage-1) + 1
1237 DO m = 1,NADJ
1238 ! Add V_i
1239 CALL WAXPY(NVAR,ONE,V(istart,m),1,Lambda(1,m),1)
1240 ! Add (H0xK_i)^T * U_i
1241 CALL HessTR_Vec ( Hes0, U(istart,m), K(istart), Tmp )
1242 CALL WAXPY(NVAR,ONE,Tmp,1,Lambda(1,m),1)
1243 END DO ! m=1:NADJ
1244 END DO
1245 ! Add H * dJac_dT_0^T * \sum(gamma_i U_i)
1246 ! Tmp holds sum gamma_i U_i
1247 IF (.NOT.Autonomous) THEN
1248 DO m = 1,NADJ
1249 Tmp(1:NVAR) = ZERO
1250 DO istage = 1, ros_S
1251 istart = NVAR*(istage-1) + 1
1252 CALL WAXPY(NVAR,ros_Gamma(istage),U(istart,m),1,Tmp,1)
1253 END DO
1254 CALL JacTR_SP_Vec(dJdT,Tmp,Tmp2)
1255 CALL WAXPY(NVAR,H,Tmp2,1,Lambda(1,m),1)
1256 END DO ! m=1:NADJ
1257 END IF ! .NOT.Autonomous
1258
1259
1260 END DO TimeLoop
1261
1262 !~~~> Save last state
1263
1264 !~~~> Succesful exit
1265 IERR = 1 !~~~> The integration was successful
1266
1267 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1268 END SUBROUTINE ros_DadjInt
1269 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1270
1271
1272
1273 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1274 SUBROUTINE ros_CadjInt ( &
1275 NADJ, Y, &
1276 Tstart, Tend, T, &
1277 AbsTol, RelTol, &
1278 !~~~> RosenbrockADJ method coefficients
1279 ros_S, ros_M, ros_E, ros_A, ros_C, &
1280 ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
1281 !~~~> Integration parameters
1282 Autonomous, VectorTol, AdjointType, &
1283 Max_no_steps, &
1284 Roundoff, Hmin, Hmax, Hstart, Hexit, &
1285 FacMin, FacMax, FacRej, FacSafe, &
1286 !~~~> Error indicator
1287 IERR )
1288 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1289 ! Template for the implementation of a generic RosenbrockADJ method
1290 ! defined by ros_S (no of stages)
1291 ! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1292 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1293
1294 IMPLICIT NONE
1295
1296 !~~~> Input: the initial condition at Tstart; Output: the solution at T
1297 INTEGER, INTENT(IN) :: NADJ
1298 KPP_REAL, INTENT(INOUT) :: Y(NVAR,NADJ)
1299 !~~~> Input: integration interval
1300 KPP_REAL, INTENT(IN) :: Tstart,Tend
1301 !~~~> Output: time at which the solution is returned (T=Tend if success)
1302 KPP_REAL, INTENT(OUT) :: T
1303 !~~~> Input: tolerances
1304 KPP_REAL, INTENT(IN) :: AbsTol(NVAR), RelTol(NVAR)
1305 !~~~> Input: The RosenbrockADJ method parameters
1306 INTEGER, INTENT(IN) :: ros_S
1307 KPP_REAL, INTENT(IN) :: ros_M(ros_S), ros_E(ros_S), &
1308 ros_Alpha(ros_S), ros_A(ros_S*(ros_S-1)/2), &
1309 ros_Gamma(ros_S), ros_C(ros_S*(ros_S-1)/2), ros_ELO
1310 LOGICAL, INTENT(IN) :: ros_NewF(ros_S)
1311 !~~~> Input: integration parameters
1312 LOGICAL, INTENT(IN) :: Autonomous, VectorTol
1313 INTEGER, INTENT(IN) :: AdjointType
1314 KPP_REAL, INTENT(IN) :: Hstart, Hmin, Hmax
1315 INTEGER, INTENT(IN) :: Max_no_steps
1316 KPP_REAL, INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
1317 !~~~> Output: last accepted step
1318 KPP_REAL, INTENT(OUT) :: Hexit
1319 !~~~> Output: Error indicator
1320 INTEGER, INTENT(OUT) :: IERR
1321 ! ~~~~ Local variables
1322 KPP_REAL :: Y0(NVAR)
1323 KPP_REAL :: Ynew(NVAR,NADJ), Fcn0(NVAR,NADJ), Fcn(NVAR,NADJ)
1324 KPP_REAL :: K(NVAR*ros_S,NADJ), dFdT(NVAR,NADJ)
1325 KPP_REAL :: Jac0(LU_NONZERO), Ghimj(LU_NONZERO)
1326 KPP_REAL :: Jac(LU_NONZERO), dJdT(LU_NONZERO)
1327 KPP_REAL :: H, Hnew, HC, HG, Fac, Tau
1328 KPP_REAL :: Err, Yerr(NVAR,NADJ)
1329 INTEGER :: Pivot(NVAR), Direction, ioffset, i, j, istage, iadj
1330 LOGICAL :: RejectLastH, RejectMoreH, Singular
1331 !~~~> Local parameters
1332 KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
1333 KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
1334 !~~~> Locally called functions
1335 ! KPP_REAL WLAMCH
1336 ! EXTERNAL WLAMCH
1337 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1338
1339
1340 !~~~> INITIAL PREPARATIONS
1341 T = Tstart
1342 Hexit = 0.0_dp
1343 H = MIN(Hstart,Hmax)
1344 IF (ABS(H) <= 10.0_dp*Roundoff) H = DeltaMin
1345
1346 IF (Tend >= Tstart) THEN
1347 Direction = +1
1348 ELSE
1349 Direction = -1
1350 END IF
1351
1352 RejectLastH=.FALSE.
1353 RejectMoreH=.FALSE.
1354
1355 OPEN(55,file='KPP_ROOT_full_cadj.dat')
1356
1357 !~~~> Time loop begins below
1358
1359 TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff <= ZERO) &
1360 .OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) )
1361
1362 IF ( Nstp > Max_no_steps ) THEN ! Too many steps
1363 CALL ros_ErrorMsg(-6,T,H,IERR)
1364 RETURN
1365 END IF
1366 IF ( ((T+0.1d0*H) == T).OR.(H <= Roundoff) ) THEN ! Step size too small
1367 CALL ros_ErrorMsg(-7,T,H,IERR)
1368 RETURN
1369 END IF
1370
1371 !~~~> Limit H if necessary to avoid going beyond Tend
1372 Hexit = H
1373 H = MIN(H,ABS(Tend-T))
1374
1375 !~~~> Interpolate forward solution
1376 CALL ros_cadj_Y( T, Y0 )
1377 !~~~> Compute the Jacobian at current time
1378 CALL JacTemplate(T, Y0, Jac0)
1379
1380 WRITE(55,55) T, H, Y0(ind_NO2), Y0(ind_O3), &
1381 Y(ind_NO2,1), Y(ind_O3,2), &
1382 Y(ind_NO2,2), Y(ind_O3,1)
1383
1384 !~~~> Compute the function derivative with respect to T
1385 IF (.NOT.Autonomous) THEN
1386 CALL ros_JacTimeDerivative ( T, Roundoff, Y0, &
1387 Jac0, dJdT )
1388 DO iadj = 1, NADJ
1389 CALL JacTR_SP_Vec(dJdT,Y(1,iadj),dFdT(1,iadj))
1390 CALL WSCAL(NVAR,(-ONE),dFdT(1,iadj),1)
1391 END DO
1392 END IF
1393
1394 !~~~> Ydot = -J^T*Y
1395 CALL WSCAL(LU_NONZERO,(-ONE),Jac0,1)
1396 DO iadj = 1, NADJ
1397 CALL JacTR_SP_Vec(Jac0,Y(1,iadj),Fcn0(1,iadj))
1398 END DO
1399
1400 !~~~> Repeat step calculation until current step accepted
1401 UntilAccepted: DO
1402
1403 CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1), &
1404 Jac0,Ghimj,Pivot,Singular)
1405 IF (Singular) THEN ! More than 5 consecutive failed decompositions
1406 CALL ros_ErrorMsg(-8,T,H,IERR)
1407 RETURN
1408 END IF
1409
1410 !~~~> Compute the stages
1411 Stage: DO istage = 1, ros_S
1412
1413 ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR)
1414 ioffset = NVAR*(istage-1)
1415
1416 ! For the 1st istage the function has been computed previously
1417 IF ( istage == 1 ) THEN
1418 DO iadj = 1, NADJ
1419 CALL WCOPY(NVAR,Fcn0(1,iadj),1,Fcn(1,iadj),1)
1420 END DO
1421 ! istage>1 and a new function evaluation is needed at the current istage
1422 ELSEIF ( ros_NewF(istage) ) THEN
1423 CALL WCOPY(NVAR*NADJ,Y,1,Ynew,1)
1424 DO j = 1, istage-1
1425 DO iadj = 1, NADJ
1426 CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), &
1427 K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1)
1428 END DO
1429 END DO
1430 Tau = T + ros_Alpha(istage)*Direction*H
1431 ! CALL FunTemplate(Tau,Ynew,Fcn)
1432 CALL ros_cadj_Y( Tau, Y0 )
1433 CALL JacTemplate(Tau, Y0, Jac)
1434 CALL WSCAL(LU_NONZERO,(-ONE),Jac,1)
1435 DO iadj = 1, NADJ
1436 CALL JacTR_SP_Vec(Jac,Ynew(1,iadj),Fcn(1,iadj))
1437 !CALL WSCAL(NVAR,(-ONE),Fcn(1,iadj),1)
1438 END DO
1439 END IF ! if istage == 1 elseif ros_NewF(istage)
1440
1441 DO iadj = 1, NADJ
1442 CALL WCOPY(NVAR,Fcn(1,iadj),1,K(ioffset+1,iadj),1)
1443 END DO
1444 DO j = 1, istage-1
1445 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H)
1446 DO iadj = 1, NADJ
1447 CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1,iadj),1, &
1448 K(ioffset+1,iadj),1)
1449 END DO
1450 END DO
1451 IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
1452 HG = Direction*H*ros_Gamma(istage)
1453 DO iadj = 1, NADJ
1454 CALL WAXPY(NVAR,HG,dFdT(1,iadj),1,K(ioffset+1,iadj),1)
1455 END DO
1456 END IF
1457 DO iadj = 1, NADJ
1458 CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1,iadj))
1459 END DO
1460
1461 END DO Stage
1462
1463
1464 !~~~> Compute the new solution
1465 DO iadj = 1, NADJ
1466 CALL WCOPY(NVAR,Y(1,iadj),1,Ynew(1,iadj),1)
1467 DO j=1,ros_S
1468 CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1)
1469 END DO
1470 END DO
1471
1472 !~~~> Compute the error estimation
1473 CALL WSCAL(NVAR*NADJ,ZERO,Yerr,1)
1474 DO j=1,ros_S
1475 DO iadj = 1, NADJ
1476 CALL WAXPY(NVAR,ros_E(j),K(NVAR*(j-1)+1,iadj),1,Yerr(1,iadj),1)
1477 END DO
1478 END DO
1479 !~~~> Max error among all adjoint components
1480 iadj = 1
1481 Err = ros_ErrorNorm ( Y(1,iadj), Ynew(1,iadj), Yerr(1,iadj), &
1482 AbsTol, RelTol, VectorTol )
1483
1484 !~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
1485 Fac = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
1486 Hnew = H*Fac
1487
1488 !~~~> Check the error magnitude and adjust step size
1489 Nstp = Nstp+1
1490 IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step
1491 Nacc = Nacc+1
1492 CALL WCOPY(NVAR*NADJ,Ynew,1,Y,1)
1493 T = T + Direction*H
1494 Hnew = MAX(Hmin,MIN(Hnew,Hmax))
1495 IF (RejectLastH) THEN ! No step size increase after a rejected step
1496 Hnew = MIN(Hnew,H)
1497 END IF
1498 RejectLastH = .FALSE.
1499 RejectMoreH = .FALSE.
1500 H = Hnew
1501 EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
1502 ELSE !~~~> Reject step
1503 IF (RejectMoreH) THEN
1504 Hnew = H*FacRej
1505 END IF
1506 RejectMoreH = RejectLastH
1507 RejectLastH = .TRUE.
1508 H = Hnew
1509 IF (Nacc >= 1) THEN
1510 Nrej = Nrej+1
1511 END IF
1512 END IF ! Err <= 1
1513
1514 END DO UntilAccepted
1515
1516 END DO TimeLoop
1517
1518 !~~~> Succesful exit
1519 IERR = 1 !~~~> The integration was successful
1520
1521 WRITE(55,55) T, H, Y0(ind_NO2), Y0(ind_O3), &
1522 Y(ind_NO2,1), Y(ind_O3,2), &
1523 Y(ind_NO2,2), Y(ind_O3,1)
1524
1525 55 FORMAT(100(E12.5,2X))
1526 CLOSE(55)
1527
1528 END SUBROUTINE ros_CadjInt
1529
1530
1531 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1532 SUBROUTINE ros_SimpleCadjInt ( &
1533 NADJ, Y, &
1534 Tstart, Tend, T, &
1535 AbsTol, RelTol, &
1536 !~~~> RosenbrockADJ method coefficients
1537 ros_S, ros_M, ros_E, ros_A, ros_C, &
1538 ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
1539 !~~~> Integration parameters
1540 Autonomous, VectorTol, AdjointType, &
1541 Max_no_steps, &
1542 Roundoff, Hmin, Hmax, Hstart, &
1543 FacMin, FacMax, FacRej, FacSafe, &
1544 !~~~> Error indicator
1545 IERR )
1546 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1547 ! Template for the implementation of a generic RosenbrockADJ method
1548 ! defined by ros_S (no of stages)
1549 ! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1550 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1551
1552 IMPLICIT NONE
1553
1554 !~~~> Input: the initial condition at Tstart; Output: the solution at T
1555 INTEGER, INTENT(IN) :: NADJ
1556 KPP_REAL, INTENT(INOUT) :: Y(NVAR,NADJ)
1557 !~~~> Input: integration interval
1558 KPP_REAL, INTENT(IN) :: Tstart,Tend
1559 !~~~> Output: time at which the solution is returned (T=Tend if success)
1560 KPP_REAL, INTENT(OUT) :: T
1561 !~~~> Input: tolerances
1562 KPP_REAL, INTENT(IN) :: AbsTol(NVAR), RelTol(NVAR)
1563 !~~~> Input: The RosenbrockADJ method parameters
1564 INTEGER, INTENT(IN) :: ros_S
1565 KPP_REAL, INTENT(IN) :: ros_M(ros_S), ros_E(ros_S), &
1566 ros_Alpha(ros_S), ros_A(ros_S*(ros_S-1)/2), &
1567 ros_Gamma(ros_S), ros_C(ros_S*(ros_S-1)/2), ros_ELO
1568 LOGICAL, INTENT(IN) :: ros_NewF(ros_S)
1569 !~~~> Input: integration parameters
1570 LOGICAL, INTENT(IN) :: Autonomous, VectorTol
1571 INTEGER, INTENT(IN) :: AdjointType
1572 KPP_REAL, INTENT(IN) :: Hstart, Hmin, Hmax
1573 INTEGER, INTENT(IN) :: Max_no_steps
1574 KPP_REAL, INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
1575 !~~~> Output: Error indicator
1576 INTEGER, INTENT(OUT) :: IERR
1577 ! ~~~~ Local variables
1578 KPP_REAL :: Y0(NVAR), Y0old(NVAR), Told
1579 KPP_REAL :: Ynew(NVAR,NADJ), Fcn0(NVAR,NADJ), Fcn(NVAR,NADJ)
1580 KPP_REAL :: K(NVAR*ros_S,NADJ), dFdT(NVAR,NADJ)
1581 KPP_REAL :: Jac0(LU_NONZERO), Ghimj(LU_NONZERO)
1582 KPP_REAL :: Jac(LU_NONZERO), dJdT(LU_NONZERO)
1583 KPP_REAL :: H, Hnew, HC, HG, Fac, Tau
1584 KPP_REAL :: Err, ghinv
1585 INTEGER :: Pivot(NVAR), Direction, ioffset, i, j, istage, iadj
1586 INTEGER :: istack
1587 LOGICAL :: RejectLastH, RejectMoreH, Singular
1588 !~~~> Local parameters
1589 KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
1590 KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
1591 !~~~> Locally called functions
1592 ! KPP_REAL WLAMCH
1593 ! EXTERNAL WLAMCH
1594 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1595
1596
1597 !~~~> INITIAL PREPARATIONS
1598
1599 IF (Tend >= Tstart) THEN
1600 Direction = -1
1601 ELSE
1602 Direction = +1
1603 END IF
1604
1605 OPEN(55,file='KPP_ROOT_smpl_cadj.dat')
1606
1607
1608 !~~~> Time loop begins below
1609 TimeLoop: DO istack = stack_ptr,2,-1
1610
1611 T = buf_T(istack)
1612 H = buf_H(istack-1)
1613 CALL WCOPY(NVAR,buf_Y(1,istack),1,Y0,1)
1614
1615 WRITE(55,55) T, H, Y0(ind_NO2), Y0(ind_O3), &
1616 Y(ind_NO2,1), Y(ind_O3,2), Y(ind_NO2,2), Y(ind_O3,1)
1617
1618 !~~~> Compute the Jacobian at current time
1619 CALL JacTemplate(T, Y0, Jac0)
1620
1621 !~~~> Compute the function derivative with respect to T
1622 IF (.NOT.Autonomous) THEN
1623 CALL ros_JacTimeDerivative ( T, Roundoff, Y0, &
1624 Jac0, dJdT )
1625 DO iadj = 1, NADJ
1626 CALL JacTR_SP_Vec(dJdT,Y(1,iadj),dFdT(1,iadj))
1627 CALL WSCAL(NVAR,(-ONE),dFdT(1,iadj),1)
1628 END DO
1629 END IF
1630
1631 !~~~> Ydot = -J^T*Y
1632 CALL WSCAL(LU_NONZERO,(-ONE),Jac0,1)
1633 DO iadj = 1, NADJ
1634 CALL JacTR_SP_Vec(Jac0,Y(1,iadj),Fcn0(1,iadj))
1635 END DO
1636
1637 !~~~> Construct Ghimj = 1/(H*ham) - Jac0
1638 CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1)
1639 CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
1640 ghinv = ONE/(Direction*H*ros_Gamma(1))
1641 DO i=1,NVAR
1642 Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv
1643 END DO
1644 !~~~> Compute LU decomposition
1645 CALL ros_Decomp( Ghimj, Pivot, j )
1646 IF (j /= 0) THEN
1647 CALL ros_ErrorMsg(-8,T,H,IERR)
1648 PRINT*,' The matrix is singular !'
1649 STOP
1650 END IF
1651
1652 !~~~> Compute the stages
1653 Stage: DO istage = 1, ros_S
1654
1655 ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR)
1656 ioffset = NVAR*(istage-1)
1657
1658 ! For the 1st istage the function has been computed previously
1659 IF ( istage == 1 ) THEN
1660 DO iadj = 1, NADJ
1661 CALL WCOPY(NVAR,Fcn0(1,iadj),1,Fcn(1,iadj),1)
1662 END DO
1663 ! istage>1 and a new function evaluation is needed at the current istage
1664 ELSEIF ( ros_NewF(istage) ) THEN
1665 CALL WCOPY(NVAR*NADJ,Y,1,Ynew,1)
1666 DO j = 1, istage-1
1667 DO iadj = 1, NADJ
1668 CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), &
1669 K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1)
1670 END DO
1671 END DO
1672 Tau = T + ros_Alpha(istage)*Direction*H
1673 CALL ros_Hermite3( buf_T(istack-1), buf_T(istack), Tau, &
1674 buf_Y(1,istack-1), buf_Y(1,istack), &
1675 buf_dY(1,istack-1), buf_dY(1,istack), Y0 )
1676 CALL JacTemplate(Tau, Y0, Jac)
1677 CALL WSCAL(LU_NONZERO,(-ONE),Jac,1)
1678 DO iadj = 1, NADJ
1679 CALL JacTR_SP_Vec(Jac,Ynew(1,iadj),Fcn(1,iadj))
1680 END DO
1681 END IF ! if istage == 1 elseif ros_NewF(istage)
1682
1683 DO iadj = 1, NADJ
1684 CALL WCOPY(NVAR,Fcn(1,iadj),1,K(ioffset+1,iadj),1)
1685 END DO
1686 DO j = 1, istage-1
1687 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H)
1688 DO iadj = 1, NADJ
1689 CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1,iadj),1, &
1690 K(ioffset+1,iadj),1)
1691 END DO
1692 END DO
1693 IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
1694 HG = Direction*H*ros_Gamma(istage)
1695 DO iadj = 1, NADJ
1696 CALL WAXPY(NVAR,HG,dFdT(1,iadj),1,K(ioffset+1,iadj),1)
1697 END DO
1698 END IF
1699 DO iadj = 1, NADJ
1700 CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1,iadj))
1701 END DO
1702
1703 END DO Stage
1704
1705
1706 !~~~> Compute the new solution
1707 DO iadj = 1, NADJ
1708 DO j=1,ros_S
1709 CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1,iadj),1,Y(1,iadj),1)
1710 END DO
1711 END DO
1712
1713 END DO TimeLoop
1714
1715 !~~~> Succesful exit
1716 IERR = 1 !~~~> The integration was successful
1717
1718 WRITE(55,55) T, H, Y0(ind_NO2), Y0(ind_O3), &
1719 Y(ind_NO2,1), Y(ind_O3,2), &
1720 Y(ind_NO2,2), Y(ind_O3,1)
1721
1722 55 FORMAT(100(E12.5,2X))
1723 CLOSE(55)
1724
1725 END SUBROUTINE ros_SimpleCadjInt
1726
1727 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1728 KPP_REAL FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
1729 AbsTol, RelTol, VectorTol )
1730 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1731 !~~~> Computes the "scaled norm" of the error vector Yerr
1732 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1733 IMPLICIT NONE
1734
1735 ! Input arguments
1736 KPP_REAL, INTENT(IN) :: Y(NVAR), Ynew(NVAR), &
1737 Yerr(NVAR), AbsTol(NVAR), RelTol(NVAR)
1738 LOGICAL, INTENT(IN) :: VectorTol
1739 ! Local variables
1740 KPP_REAL :: Err, Scale, Ymax
1741 INTEGER :: i
1742 KPP_REAL, PARAMETER :: ZERO = 0.0d0
1743
1744 Err = ZERO
1745 DO i=1,NVAR
1746 Ymax = MAX(ABS(Y(i)),ABS(Ynew(i)))
1747 IF (VectorTol) THEN
1748 Scale = AbsTol(i)+RelTol(i)*Ymax
1749 ELSE
1750 Scale = AbsTol(1)+RelTol(1)*Ymax
1751 END IF
1752 Err = Err+(Yerr(i)/Scale)**2
1753 END DO
1754 Err = SQRT(Err/NVAR)
1755
1756 ros_ErrorNorm = Err
1757
1758 END FUNCTION ros_ErrorNorm
1759
1760
1761 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1762 SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT )
1763 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1764 !~~~> The time partial derivative of the function by finite differences
1765 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1766 IMPLICIT NONE
1767
1768 !~~~> Input arguments
1769 KPP_REAL, INTENT(IN) :: T, Roundoff, Y(NVAR), Fcn0(NVAR)
1770 !~~~> Output arguments
1771 KPP_REAL, INTENT(OUT) :: dFdT(NVAR)
1772 !~~~> Local variables
1773 KPP_REAL :: Delta
1774 KPP_REAL, PARAMETER :: ONE = 1.0d0, DeltaMin = 1.0d-6
1775
1776 Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
1777 CALL FunTemplate(T+Delta,Y,dFdT)
1778 CALL WAXPY(NVAR,(-ONE),Fcn0,1,dFdT,1)
1779 CALL WSCAL(NVAR,(ONE/Delta),dFdT,1)
1780
1781 END SUBROUTINE ros_FunTimeDerivative
1782
1783
1784 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1785 SUBROUTINE ros_JacTimeDerivative ( T, Roundoff, Y, &
1786 Jac0, dJdT )
1787 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1788 !~~~> The time partial derivative of the Jacobian by finite differences
1789 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1790 IMPLICIT NONE
1791
1792 !~~~> Input arguments
1793 KPP_REAL, INTENT(IN) :: T, Roundoff, Y(NVAR), Jac0(LU_NONZERO)
1794 !~~~> Output arguments
1795 KPP_REAL, INTENT(OUT) :: dJdT(LU_NONZERO)
1796 !~~~> Local variables
1797 KPP_REAL Delta
1798 KPP_REAL, PARAMETER :: ONE = 1.0d0, DeltaMin = 1.0d-6
1799
1800 Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
1801 CALL JacTemplate(T+Delta,Y,dJdT)
1802 CALL WAXPY(LU_NONZERO,(-ONE),Jac0,1,dJdT,1)
1803 CALL WSCAL(LU_NONZERO,(ONE/Delta),dJdT,1)
1804
1805 END SUBROUTINE ros_JacTimeDerivative
1806
1807
1808 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1809 SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, &
1810 Jac0, Ghimj, Pivot, Singular )
1811 ! --- --- --- --- --- --- --- --- --- --- --- --- ---
1812 ! Prepares the LHS matrix for stage calculations
1813 ! 1. Construct Ghimj = 1/(H*ham) - Jac0
1814 ! "(Gamma H) Inverse Minus Jacobian"
1815 ! 2. Repeat LU decomposition of Ghimj until successful.
1816 ! -half the step size if LU decomposition fails and retry
1817 ! -exit after 5 consecutive fails
1818 ! --- --- --- --- --- --- --- --- --- --- --- --- ---
1819 IMPLICIT NONE
1820
1821 !~~~> Input arguments
1822 KPP_REAL, INTENT(IN) :: gam, Jac0(LU_NONZERO)
1823 INTEGER, INTENT(IN) :: Direction
1824 !~~~> Output arguments
1825 KPP_REAL, INTENT(OUT) :: Ghimj(LU_NONZERO)
1826 LOGICAL, INTENT(OUT) :: Singular
1827 INTEGER, INTENT(OUT) :: Pivot(NVAR)
1828 !~~~> Inout arguments
1829 KPP_REAL, INTENT(INOUT) :: H ! step size is decreased when LU fails
1830 !~~~> Local variables
1831 INTEGER :: i, ising, Nconsecutive
1832 KPP_REAL :: ghinv
1833 KPP_REAL, PARAMETER :: ONE = 1.0d0, HALF = 0.5d0
1834
1835 Nconsecutive = 0
1836 Singular = .TRUE.
1837
1838 DO WHILE (Singular)
1839
1840 !~~~> Construct Ghimj = 1/(H*ham) - Jac0
1841 CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1)
1842 CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
1843 ghinv = ONE/(Direction*H*gam)
1844 DO i=1,NVAR
1845 Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv
1846 END DO
1847 !~~~> Compute LU decomposition
1848 CALL ros_Decomp( Ghimj, Pivot, ising )
1849 IF (ising == 0) THEN
1850 !~~~> If successful done
1851 Singular = .FALSE.
1852 ELSE ! ising .ne. 0
1853 !~~~> If unsuccessful half the step size; if 5 consecutive fails then return
1854 Nsng = Nsng+1
1855 Nconsecutive = Nconsecutive+1
1856 Singular = .TRUE.
1857 PRINT*,'Warning: LU Decomposition returned ising = ',ising
1858 IF (Nconsecutive <= 5) THEN ! Less than 5 consecutive failed decompositions
1859 H = H*HALF
1860 ELSE ! More than 5 consecutive failed decompositions
1861 RETURN
1862 END IF ! Nconsecutive
1863 END IF ! ising
1864
1865 END DO ! WHILE Singular
1866
1867 END SUBROUTINE ros_PrepareMatrix
1868
1869
1870 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1871 SUBROUTINE ros_Decomp( A, Pivot, ising )
1872 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1873 ! Template for the LU decomposition
1874 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1875 IMPLICIT NONE
1876 !~~~> Inout variables
1877 KPP_REAL, INTENT(INOUT) :: A(LU_NONZERO)
1878 !~~~> Output variables
1879 INTEGER, INTENT(OUT) :: Pivot(NVAR), ising
1880
1881 CALL KppDecomp ( A, ising )
1882 !~~~> Note: for a full matrix use Lapack:
1883 ! CALL DGETRF( NVAR, NVAR, A, NVAR, Pivot, ising )
1884 Pivot(1) = 1
1885
1886 Ndec = Ndec + 1
1887
1888 END SUBROUTINE ros_Decomp
1889
1890
1891 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1892 SUBROUTINE ros_Solve( C, A, Pivot, b )
1893 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1894 ! Template for the forward/backward substitution (using pre-computed LU decomposition)
1895 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1896 IMPLICIT NONE
1897 !~~~> Input variables
1898 CHARACTER, INTENT(IN) :: C
1899 KPP_REAL, INTENT(IN) :: A(LU_NONZERO)
1900 INTEGER, INTENT(IN) :: Pivot(NVAR)
1901 !~~~> InOut variables
1902 KPP_REAL, INTENT(INOUT) :: b(NVAR)
1903
1904 SELECT CASE (C)
1905 CASE ('N')
1906 CALL KppSolve( A, b )
1907 CASE ('T')
1908 CALL KppSolveTR( A, b, b )
1909 CASE DEFAULT
1910 PRINT*,'Unknown C = (',C,') in ros_Solve'
1911 STOP
1912 END SELECT
1913 !~~~> Note: for a full matrix use Lapack:
1914 ! NRHS = 1
1915 ! CALL DGETRS( C, NVAR , NRHS, A, NVAR, Pivot, b, NVAR, INFO )
1916
1917 Nsol = Nsol+1
1918
1919 END SUBROUTINE ros_Solve
1920
1921
1922 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1923 SUBROUTINE ros_cadj_Y( T, Y )
1924 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1925 ! Finds the solution Y at T by interpolating the stored forward trajectory
1926 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1927 IMPLICIT NONE
1928 !~~~> Input variables
1929 KPP_REAL, INTENT(IN) :: T
1930 !~~~> Output variables
1931 KPP_REAL, INTENT(OUT) :: Y(NVAR)
1932 !~~~> Local variables
1933 INTEGER :: i, j
1934 KPP_REAL, PARAMETER :: ONE = 1.0d0
1935
1936 ! buf_H, buf_T, buf_Y, buf_dY, buf_d2Y
1937
1938 IF( (T < buf_T(1)).OR.(T> buf_T(stack_ptr)) ) THEN
1939 PRINT*,'Cannot locate solution at T = ',T
1940 PRINT*,'Stored trajectory is between Tstart = ',buf_T(1)
1941 PRINT*,' and Tend = ',buf_T(stack_ptr)
1942 STOP
1943 END IF
1944 DO i = 1, stack_ptr-1
1945 IF( (T>= buf_T(i)).AND.(T<= buf_T(i+1)) ) EXIT
1946 END DO
1947
1948
1949 IF (.FALSE.) THEN
1950
1951 CALL ros_Hermite5( buf_T(i), buf_T(i+1), T, &
1952 buf_Y(1,i), buf_Y(1,i+1), &
1953 buf_dY(1,i), buf_dY(1,i+1), &
1954 buf_d2Y(1,i), buf_d2Y(1,i+1), Y )
1955
1956 ELSE
1957
1958 CALL ros_Hermite3( buf_T(i), buf_T(i+1), T, &
1959 buf_Y(1,i), buf_Y(1,i+1), &
1960 buf_dY(1,i), buf_dY(1,i+1), &
1961 Y )
1962
1963
1964 END IF
1965
1966 END SUBROUTINE ros_cadj_Y
1967
1968
1969 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1970 SUBROUTINE ros_Hermite3( a, b, T, Ya, Yb, Ja, Jb, Y )
1971 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1972 ! Template for Hermite interpolation of order 5 on the interval [a,b]
1973 ! P = c(1) + c(2)*(x-a) + ... + c(4)*(x-a)^3
1974 ! P[a,b] = [Ya,Yb], P'[a,b] = [Ja,Jb]
1975 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1976 IMPLICIT NONE
1977 !~~~> Input variables
1978 KPP_REAL, INTENT(IN) :: a, b, T, Ya(NVAR), Yb(NVAR)
1979 KPP_REAL, INTENT(IN) :: Ja(NVAR), Jb(NVAR)
1980 !~~~> Output variables
1981 KPP_REAL, INTENT(OUT) :: Y(NVAR)
1982 !~~~> Local variables
1983 KPP_REAL :: Tau, amb(3), C(NVAR,4)
1984 KPP_REAL, PARAMETER :: ZERO = 0.0d0
1985 INTEGER :: i, j
1986
1987 amb(1) = 1.0d0/(a-b)
1988 DO i=2,3
1989 amb(i) = amb(i-1)*amb(1)
1990 END DO
1991
1992
1993 ! c(1) = ya;
1994 CALL WCOPY(NVAR,Ya,1,C(1,1),1)
1995 ! c(2) = ja;
1996 CALL WCOPY(NVAR,Ja,1,C(1,2),1)
1997 ! c(3) = 2/(a-b)*ja + 1/(a-b)*jb - 3/(a - b)^2*ya + 3/(a - b)^2*yb ;
1998 CALL WCOPY(NVAR,Ya,1,C(1,3),1)
1999 CALL WSCAL(NVAR,-3.0*amb(2),C(1,3),1)
2000 CALL WAXPY(NVAR,3.0*amb(2),Yb,1,C(1,3),1)
2001 CALL WAXPY(NVAR,2.0*amb(1),Ja,1,C(1,3),1)
2002 CALL WAXPY(NVAR,amb(1),Jb,1,C(1,3),1)
2003 ! c(4) = 1/(a-b)^2*ja + 1/(a-b)^2*jb - 2/(a-b)^3*ya + 2/(a-b)^3*yb ;
2004 CALL WCOPY(NVAR,Ya,1,C(1,4),1)
2005 CALL WSCAL(NVAR,-2.0*amb(3),C(1,4),1)
2006 CALL WAXPY(NVAR,2.0*amb(3),Yb,1,C(1,4),1)
2007 CALL WAXPY(NVAR,amb(2),Ja,1,C(1,4),1)
2008 CALL WAXPY(NVAR,amb(2),Jb,1,C(1,4),1)
2009
2010 Tau = T - a
2011 CALL WCOPY(NVAR,C(1,4),1,Y,1)
2012 CALL WSCAL(NVAR,Tau**3,Y,1)
2013 DO j = 3,1,-1
2014 CALL WAXPY(NVAR,TAU**(j-1),C(1,j),1,Y,1)
2015 END DO
2016
2017 END SUBROUTINE ros_Hermite3
2018
2019 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2020 SUBROUTINE ros_Hermite5( a, b, T, Ya, Yb, Ja, Jb, Ha, Hb, Y )
2021 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2022 ! Template for Hermite interpolation of order 5 on the interval [a,b]
2023 ! P = c(1) + c(2)*(x-a) + ... + c(6)*(x-a)^5
2024 ! P[a,b] = [Ya,Yb], P'[a,b] = [Ja,Jb], P"[a,b] = [Ha,Hb]
2025 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2026 IMPLICIT NONE
2027 !~~~> Input variables
2028 KPP_REAL, INTENT(IN) :: a, b, T, Ya(NVAR), Yb(NVAR)
2029 KPP_REAL, INTENT(IN) :: Ja(NVAR), Jb(NVAR), Ha(NVAR), Hb(NVAR)
2030 !~~~> Output variables
2031 KPP_REAL, INTENT(OUT) :: Y(NVAR)
2032 !~~~> Local variables
2033 KPP_REAL :: Tau, amb(5), C(NVAR,6)
2034 KPP_REAL, PARAMETER :: ZERO = 0.0d0, HALF = 0.5d0
2035 INTEGER :: i, j
2036
2037 amb(1) = 1.0d0/(a-b)
2038 DO i=2,5
2039 amb(i) = amb(i-1)*amb(1)
2040 END DO
2041
2042 ! c(1) = ya;
2043 CALL WCOPY(NVAR,Ya,1,C(1,1),1)
2044 ! c(2) = ja;
2045 CALL WCOPY(NVAR,Ja,1,C(1,2),1)
2046 ! c(3) = ha/2;
2047 CALL WCOPY(NVAR,Ha,1,C(1,3),1)
2048 CALL WSCAL(NVAR,HALF,C(1,3),1)
2049
2050 ! c(4) = 10*amb(3)*ya - 10*amb(3)*yb - 6*amb(2)*ja - 4*amb(2)*jb + 1.5*amb(1)*ha - 0.5*amb(1)*hb ;
2051 CALL WCOPY(NVAR,Ya,1,C(1,4),1)
2052 CALL WSCAL(NVAR,10.0*amb(3),C(1,4),1)
2053 CALL WAXPY(NVAR,-10.0*amb(3),Yb,1,C(1,4),1)
2054 CALL WAXPY(NVAR,-6.0*amb(2),Ja,1,C(1,4),1)
2055 CALL WAXPY(NVAR,-4.0*amb(2),Jb,1,C(1,4),1)
2056 CALL WAXPY(NVAR, 1.5*amb(1),Ha,1,C(1,4),1)
2057 CALL WAXPY(NVAR,-0.5*amb(1),Hb,1,C(1,4),1)
2058
2059 ! c(5) = 15*amb(4)*ya - 15*amb(4)*yb - 8.*amb(3)*ja - 7*amb(3)*jb + 1.5*amb(2)*ha - 1*amb(2)*hb ;
2060 CALL WCOPY(NVAR,Ya,1,C(1,5),1)
2061 CALL WSCAL(NVAR, 15.0*amb(4),C(1,5),1)
2062 CALL WAXPY(NVAR,-15.0*amb(4),Yb,1,C(1,5),1)
2063 CALL WAXPY(NVAR,-8.0*amb(3),Ja,1,C(1,5),1)
2064 CALL WAXPY(NVAR,-7.0*amb(3),Jb,1,C(1,5),1)
2065 CALL WAXPY(NVAR,1.5*amb(2),Ha,1,C(1,5),1)
2066 CALL WAXPY(NVAR,-amb(2),Hb,1,C(1,5),1)
2067
2068 ! c(6) = 6*amb(5)*ya - 6*amb(5)*yb - 3.*amb(4)*ja - 3.*amb(4)*jb + 0.5*amb(3)*ha -0.5*amb(3)*hb ;
2069 CALL WCOPY(NVAR,Ya,1,C(1,6),1)
2070 CALL WSCAL(NVAR, 6.0*amb(5),C(1,6),1)
2071 CALL WAXPY(NVAR,-6.0*amb(5),Yb,1,C(1,6),1)
2072 CALL WAXPY(NVAR,-3.0*amb(4),Ja,1,C(1,6),1)
2073 CALL WAXPY(NVAR,-3.0*amb(4),Jb,1,C(1,6),1)
2074 CALL WAXPY(NVAR, 0.5*amb(3),Ha,1,C(1,6),1)
2075 CALL WAXPY(NVAR,-0.5*amb(3),Hb,1,C(1,6),1)
2076
2077 Tau = T - a
2078 CALL WCOPY(NVAR,C(1,6),1,Y,1)
2079 DO j = 5,1,-1
2080 CALL WSCAL(NVAR,Tau,Y,1)
2081 CALL WAXPY(NVAR,ONE,C(1,j),1,Y,1)
2082 END DO
2083
2084 END SUBROUTINE ros_Hermite5
2085
2086 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2087 SUBROUTINE Ros2 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
2088 ros_Gamma,ros_NewF,ros_ELO,ros_Name)
2089 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2090 ! --- AN L-STABLE METHOD, 2 stages, order 2
2091 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2092
2093 IMPLICIT NONE
2094
2095 INTEGER, PARAMETER :: S = 2
2096 INTEGER, INTENT(OUT) :: ros_S
2097 KPP_REAL, DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
2098 KPP_REAL, DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
2099 KPP_REAL, INTENT(OUT) :: ros_ELO
2100 LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
2101 CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
2102 DOUBLE PRECISION g
2103
2104 g = 1.0d0 + 1.0d0/SQRT(2.0d0)
2105
2106 !~~~> Name of the method
2107 ros_Name = 'ROS-2'
2108 !~~~> Number of stages
2109 ros_S = S
2110
2111 !~~~> The coefficient matrices A and C are strictly lower triangular.
2112 ! The lower triangular (subdiagonal) elements are stored in row-wise order:
2113 ! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2114 ! The general mapping formula is:
2115 ! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2116 ! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2117
2118 ros_A(1) = (1.d0)/g
2119 ros_C(1) = (-2.d0)/g
2120 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2121 ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2122 ros_NewF(1) = .TRUE.
2123 ros_NewF(2) = .TRUE.
2124 !~~~> M_i = Coefficients for new step solution
2125 ros_M(1)= (3.d0)/(2.d0*g)
2126 ros_M(2)= (1.d0)/(2.d0*g)
2127 ! E_i = Coefficients for error estimator
2128 ros_E(1) = 1.d0/(2.d0*g)
2129 ros_E(2) = 1.d0/(2.d0*g)
2130 !~~~> ros_ELO = estimator of local order - the minimum between the
2131 ! main and the embedded scheme orders plus one
2132 ros_ELO = 2.0d0
2133 !~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2134 ros_Alpha(1) = 0.0d0
2135 ros_Alpha(2) = 1.0d0
2136 !~~~> Gamma_i = \sum_j gamma_{i,j}
2137 ros_Gamma(1) = g
2138 ros_Gamma(2) =-g
2139
2140 END SUBROUTINE Ros2
2141
2142
2143 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2144 SUBROUTINE Ros3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
2145 ros_Gamma,ros_NewF,ros_ELO,ros_Name)
2146 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2147 ! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
2148 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2149
2150 IMPLICIT NONE
2151
2152 INTEGER, PARAMETER :: S = 3
2153 INTEGER, INTENT(OUT) :: ros_S
2154 KPP_REAL, DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
2155 KPP_REAL, DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
2156 KPP_REAL, INTENT(OUT) :: ros_ELO
2157 LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
2158 CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
2159
2160 !~~~> Name of the method
2161 ros_Name = 'ROS-3'
2162 !~~~> Number of stages
2163 ros_S = S
2164
2165 !~~~> The coefficient matrices A and C are strictly lower triangular.
2166 ! The lower triangular (subdiagonal) elements are stored in row-wise order:
2167 ! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2168 ! The general mapping formula is:
2169 ! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2170 ! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2171
2172 ros_A(1)= 1.d0
2173 ros_A(2)= 1.d0
2174 ros_A(3)= 0.d0
2175
2176 ros_C(1) = -0.10156171083877702091975600115545d+01
2177 ros_C(2) = 0.40759956452537699824805835358067d+01
2178 ros_C(3) = 0.92076794298330791242156818474003d+01
2179 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2180 ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2181 ros_NewF(1) = .TRUE.
2182 ros_NewF(2) = .TRUE.
2183 ros_NewF(3) = .FALSE.
2184 !~~~> M_i = Coefficients for new step solution
2185 ros_M(1) = 0.1d+01
2186 ros_M(2) = 0.61697947043828245592553615689730d+01
2187 ros_M(3) = -0.42772256543218573326238373806514d+00
2188 ! E_i = Coefficients for error estimator
2189 ros_E(1) = 0.5d+00
2190 ros_E(2) = -0.29079558716805469821718236208017d+01
2191 ros_E(3) = 0.22354069897811569627360909276199d+00
2192 !~~~> ros_ELO = estimator of local order - the minimum between the
2193 ! main and the embedded scheme orders plus 1
2194 ros_ELO = 3.0d0
2195 !~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2196 ros_Alpha(1)= 0.0d+00
2197 ros_Alpha(2)= 0.43586652150845899941601945119356d+00
2198 ros_Alpha(3)= 0.43586652150845899941601945119356d+00
2199 !~~~> Gamma_i = \sum_j gamma_{i,j}
2200 ros_Gamma(1)= 0.43586652150845899941601945119356d+00
2201 ros_Gamma(2)= 0.24291996454816804366592249683314d+00
2202 ros_Gamma(3)= 0.21851380027664058511513169485832d+01
2203
2204 END SUBROUTINE Ros3
2205
2206 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2207
2208
2209 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2210 SUBROUTINE Ros4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
2211 ros_Gamma,ros_NewF,ros_ELO,ros_Name)
2212 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2213 ! L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
2214 ! L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
2215 !
2216 ! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
2217 ! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2218 ! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2219 ! SPRINGER-VERLAG (1990)
2220 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2221
2222 IMPLICIT NONE
2223
2224 INTEGER, PARAMETER :: S=4
2225 INTEGER, INTENT(OUT) :: ros_S
2226 KPP_REAL, DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
2227 KPP_REAL, DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
2228 KPP_REAL, INTENT(OUT) :: ros_ELO
2229 LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
2230 CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
2231 DOUBLE PRECISION g
2232
2233 !~~~> Name of the method
2234 ros_Name = 'ROS-4'
2235 !~~~> Number of stages
2236 ros_S = S
2237
2238 !~~~> The coefficient matrices A and C are strictly lower triangular.
2239 ! The lower triangular (subdiagonal) elements are stored in row-wise order:
2240 ! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2241 ! The general mapping formula is:
2242 ! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2243 ! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2244
2245 ros_A(1) = 0.2000000000000000d+01
2246 ros_A(2) = 0.1867943637803922d+01
2247 ros_A(3) = 0.2344449711399156d+00
2248 ros_A(4) = ros_A(2)
2249 ros_A(5) = ros_A(3)
2250 ros_A(6) = 0.0D0
2251
2252 ros_C(1) =-0.7137615036412310d+01
2253 ros_C(2) = 0.2580708087951457d+01
2254 ros_C(3) = 0.6515950076447975d+00
2255 ros_C(4) =-0.2137148994382534d+01
2256 ros_C(5) =-0.3214669691237626d+00
2257 ros_C(6) =-0.6949742501781779d+00
2258 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2259 ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2260 ros_NewF(1) = .TRUE.
2261 ros_NewF(2) = .TRUE.
2262 ros_NewF(3) = .TRUE.
2263 ros_NewF(4) = .FALSE.
2264 !~~~> M_i = Coefficients for new step solution
2265 ros_M(1) = 0.2255570073418735d+01
2266 ros_M(2) = 0.2870493262186792d+00
2267 ros_M(3) = 0.4353179431840180d+00
2268 ros_M(4) = 0.1093502252409163d+01
2269 !~~~> E_i = Coefficients for error estimator
2270 ros_E(1) =-0.2815431932141155d+00
2271 ros_E(2) =-0.7276199124938920d-01
2272 ros_E(3) =-0.1082196201495311d+00
2273 ros_E(4) =-0.1093502252409163d+01
2274 !~~~> ros_ELO = estimator of local order - the minimum between the
2275 ! main and the embedded scheme orders plus 1
2276 ros_ELO = 4.0d0
2277 !~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2278 ros_Alpha(1) = 0.D0
2279 ros_Alpha(2) = 0.1145640000000000d+01
2280 ros_Alpha(3) = 0.6552168638155900d+00
2281 ros_Alpha(4) = ros_Alpha(3)
2282 !~~~> Gamma_i = \sum_j gamma_{i,j}
2283 ros_Gamma(1) = 0.5728200000000000d+00
2284 ros_Gamma(2) =-0.1769193891319233d+01
2285 ros_Gamma(3) = 0.7592633437920482d+00
2286 ros_Gamma(4) =-0.1049021087100450d+00
2287
2288 END SUBROUTINE Ros4
2289
2290 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2291 SUBROUTINE Rodas3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
2292 ros_Gamma,ros_NewF,ros_ELO,ros_Name)
2293 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2294 ! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
2295 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2296
2297 IMPLICIT NONE
2298
2299 INTEGER, PARAMETER :: S=4
2300 INTEGER, INTENT(OUT) :: ros_S
2301 KPP_REAL, DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
2302 KPP_REAL, DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
2303 KPP_REAL, INTENT(OUT) :: ros_ELO
2304 LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
2305 CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
2306 DOUBLE PRECISION g
2307
2308 !~~~> Name of the method
2309 ros_Name = 'RODAS-3'
2310 !~~~> Number of stages
2311 ros_S = S
2312
2313 !~~~> The coefficient matrices A and C are strictly lower triangular.
2314 ! The lower triangular (subdiagonal) elements are stored in row-wise order:
2315 ! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2316 ! The general mapping formula is:
2317 ! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2318 ! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2319
2320 ros_A(1) = 0.0d+00
2321 ros_A(2) = 2.0d+00
2322 ros_A(3) = 0.0d+00
2323 ros_A(4) = 2.0d+00
2324 ros_A(5) = 0.0d+00
2325 ros_A(6) = 1.0d+00
2326
2327 ros_C(1) = 4.0d+00
2328 ros_C(2) = 1.0d+00
2329 ros_C(3) =-1.0d+00
2330 ros_C(4) = 1.0d+00
2331 ros_C(5) =-1.0d+00
2332 ros_C(6) =-(8.0d+00/3.0d+00)
2333
2334 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2335 ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2336 ros_NewF(1) = .TRUE.
2337 ros_NewF(2) = .FALSE.
2338 ros_NewF(3) = .TRUE.
2339 ros_NewF(4) = .TRUE.
2340 !~~~> M_i = Coefficients for new step solution
2341 ros_M(1) = 2.0d+00
2342 ros_M(2) = 0.0d+00
2343 ros_M(3) = 1.0d+00
2344 ros_M(4) = 1.0d+00
2345 !~~~> E_i = Coefficients for error estimator
2346 ros_E(1) = 0.0d+00
2347 ros_E(2) = 0.0d+00
2348 ros_E(3) = 0.0d+00
2349 ros_E(4) = 1.0d+00
2350 !~~~> ros_ELO = estimator of local order - the minimum between the
2351 ! main and the embedded scheme orders plus 1
2352 ros_ELO = 3.0d+00
2353 !~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2354 ros_Alpha(1) = 0.0d+00
2355 ros_Alpha(2) = 0.0d+00
2356 ros_Alpha(3) = 1.0d+00
2357 ros_Alpha(4) = 1.0d+00
2358 !~~~> Gamma_i = \sum_j gamma_{i,j}
2359 ros_Gamma(1) = 0.5d+00
2360 ros_Gamma(2) = 1.5d+00
2361 ros_Gamma(3) = 0.0d+00
2362 ros_Gamma(4) = 0.0d+00
2363
2364 END SUBROUTINE Rodas3
2365
2366 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2367 SUBROUTINE Rodas4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
2368 ros_Gamma,ros_NewF,ros_ELO,ros_Name)
2369 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2370 ! STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
2371 !
2372 ! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
2373 ! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2374 ! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2375 ! SPRINGER-VERLAG (1996)
2376 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2377
2378 IMPLICIT NONE
2379
2380 INTEGER, PARAMETER :: S=6
2381 INTEGER, INTENT(OUT) :: ros_S
2382 KPP_REAL, DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
2383 KPP_REAL, DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
2384 KPP_REAL, INTENT(OUT) :: ros_ELO
2385 LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
2386 CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
2387 DOUBLE PRECISION g
2388
2389 !~~~> Name of the method
2390 ros_Name = 'RODAS-4'
2391 !~~~> Number of stages
2392 ros_S = S
2393
2394 !~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2395 ros_Alpha(1) = 0.000d0
2396 ros_Alpha(2) = 0.386d0
2397 ros_Alpha(3) = 0.210d0
2398 ros_Alpha(4) = 0.630d0
2399 ros_Alpha(5) = 1.000d0
2400 ros_Alpha(6) = 1.000d0
2401
2402 !~~~> Gamma_i = \sum_j gamma_{i,j}
2403 ros_Gamma(1) = 0.2500000000000000d+00
2404 ros_Gamma(2) =-0.1043000000000000d+00
2405 ros_Gamma(3) = 0.1035000000000000d+00
2406 ros_Gamma(4) =-0.3620000000000023d-01
2407 ros_Gamma(5) = 0.0d0
2408 ros_Gamma(6) = 0.0d0
2409
2410 !~~~> The coefficient matrices A and C are strictly lower triangular.
2411 ! The lower triangular (subdiagonal) elements are stored in row-wise order:
2412 ! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2413 ! The general mapping formula is: A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2414 ! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2415
2416 ros_A(1) = 0.1544000000000000d+01
2417 ros_A(2) = 0.9466785280815826d+00
2418 ros_A(3) = 0.2557011698983284d+00
2419 ros_A(4) = 0.3314825187068521d+01
2420 ros_A(5) = 0.2896124015972201d+01
2421 ros_A(6) = 0.9986419139977817d+00
2422 ros_A(7) = 0.1221224509226641d+01
2423 ros_A(8) = 0.6019134481288629d+01
2424 ros_A(9) = 0.1253708332932087d+02
2425 ros_A(10) =-0.6878860361058950d+00
2426 ros_A(11) = ros_A(7)
2427 ros_A(12) = ros_A(8)
2428 ros_A(13) = ros_A(9)
2429 ros_A(14) = ros_A(10)
2430 ros_A(15) = 1.0d+00
2431
2432 ros_C(1) =-0.5668800000000000d+01
2433 ros_C(2) =-0.2430093356833875d+01
2434 ros_C(3) =-0.2063599157091915d+00
2435 ros_C(4) =-0.1073529058151375d+00
2436 ros_C(5) =-0.9594562251023355d+01
2437 ros_C(6) =-0.2047028614809616d+02
2438 ros_C(7) = 0.7496443313967647d+01
2439 ros_C(8) =-0.1024680431464352d+02
2440 ros_C(9) =-0.3399990352819905d+02
2441 ros_C(10) = 0.1170890893206160d+02
2442 ros_C(11) = 0.8083246795921522d+01
2443 ros_C(12) =-0.7981132988064893d+01
2444 ros_C(13) =-0.3152159432874371d+02
2445 ros_C(14) = 0.1631930543123136d+02
2446 ros_C(15) =-0.6058818238834054d+01
2447
2448 !~~~> M_i = Coefficients for new step solution
2449 ros_M(1) = ros_A(7)
2450 ros_M(2) = ros_A(8)
2451 ros_M(3) = ros_A(9)
2452 ros_M(4) = ros_A(10)
2453 ros_M(5) = 1.0d+00
2454 ros_M(6) = 1.0d+00
2455
2456 !~~~> E_i = Coefficients for error estimator
2457 ros_E(1) = 0.0d+00
2458 ros_E(2) = 0.0d+00
2459 ros_E(3) = 0.0d+00
2460 ros_E(4) = 0.0d+00
2461 ros_E(5) = 0.0d+00
2462 ros_E(6) = 1.0d+00
2463
2464 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2465 ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2466 ros_NewF(1) = .TRUE.
2467 ros_NewF(2) = .TRUE.
2468 ros_NewF(3) = .TRUE.
2469 ros_NewF(4) = .TRUE.
2470 ros_NewF(5) = .TRUE.
2471 ros_NewF(6) = .TRUE.
2472
2473 !~~~> ros_ELO = estimator of local order - the minimum between the
2474 ! main and the embedded scheme orders plus 1
2475 ros_ELO = 4.0d0
2476
2477 END SUBROUTINE Rodas4
2478
2479
2480
2481 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2482 SUBROUTINE FunTemplate( T, Y, Ydot )
2483 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2484 ! Template for the ODE function call.
2485 ! Updates the rate coefficients (and possibly the fixed species) at each call
2486 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2487
2488 !~~~> Input variables
2489 KPP_REAL T, Y(NVAR)
2490 !~~~> Output variables
2491 KPP_REAL Ydot(NVAR)
2492 !~~~> Local variables
2493 KPP_REAL Told
2494
2495 Told = TIME
2496 TIME = T
2497 CALL Update_SUN()
2498 CALL Update_RCONST()
2499 CALL Fun( Y, FIX, RCONST, Ydot )
2500 TIME = Told
2501
2502 Nfun = Nfun+1
2503
2504 END SUBROUTINE FunTemplate
2505
2506
2507 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2508 SUBROUTINE JacTemplate( T, Y, Jcb )
2509 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2510 ! Template for the ODE Jacobian call.
2511 ! Updates the rate coefficients (and possibly the fixed species) at each call
2512 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2513
2514 !~~~> Input variables
2515 KPP_REAL T, Y(NVAR)
2516 !~~~> Output variables
2517 KPP_REAL Jcb(LU_NONZERO)
2518 !~~~> Local variables
2519 KPP_REAL Told
2520
2521 Told = TIME
2522 TIME = T
2523 CALL Update_SUN()
2524 CALL Update_RCONST()
2525 CALL Jac_SP( Y, FIX, RCONST, Jcb )
2526 TIME = Told
2527
2528 Njac = Njac+1
2529
2530 END SUBROUTINE JacTemplate
2531
2532
2533 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2534 SUBROUTINE HessTemplate( T, Y, Hes )
2535 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2536 ! Template for the ODE Hessian call.
2537 ! Updates the rate coefficients (and possibly the fixed species) at each call
2538 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2539
2540 !~~~> Input variables
2541 KPP_REAL T, Y(NVAR)
2542 !~~~> Output variables
2543 KPP_REAL Hes(NHESS)
2544 !~~~> Local variables
2545 KPP_REAL Told
2546
2547 Told = TIME
2548 TIME = T
2549 CALL Update_SUN()
2550 CALL Update_RCONST()
2551 CALL Hessian( Y, FIX, RCONST, Hes )
2552 TIME = Told
2553
2554 END SUBROUTINE HessTemplate
2555
2556 END MODULE KPP_ROOT_Integrator
2557
Mirror of the official repository at https://github.com/wrf-model/WRF