1 | SUBROUTINE ros2_cts_adj(N,T,Tnext,Hmin,Hmax,Hstart, |
---|
2 | + y,Lambda,Fix,Rconst,AbsTol,RelTol, |
---|
3 | + Info) |
---|
4 | INCLUDE 'KPP_ROOT_params.h' |
---|
5 | INCLUDE 'KPP_ROOT_global.h' |
---|
6 | |
---|
7 | C INPUT ARGUMENTS: |
---|
8 | C y = Vector of (NVAR) concentrations, contains the |
---|
9 | C initial values on input |
---|
10 | C [T, Tnext] = the integration interval |
---|
11 | C Hmin, Hmax = lower and upper bounds for the selected step-size. |
---|
12 | C Note that for Step = Hmin the current computed |
---|
13 | C solution is unconditionally accepted by the error |
---|
14 | C control mechanism. |
---|
15 | C AbsTol, RelTol = (NVAR) dimensional vectors of |
---|
16 | C componentwise absolute and relative tolerances. |
---|
17 | C FUN = name of routine of derivatives. KPP syntax. |
---|
18 | C See the header below. |
---|
19 | C JAC_SP = name of routine that computes the Jacobian, in |
---|
20 | C sparse format. KPP syntax. See the header below. |
---|
21 | C Info(1) = 1 for autonomous system |
---|
22 | C = 0 for nonautonomous system |
---|
23 | C Info(2) = 1 for third order embedded formula |
---|
24 | C = 0 for first order embedded formula |
---|
25 | C |
---|
26 | C Note: Stage 3 used to build strongly A-stable order 3 formula for error control |
---|
27 | C Embed3 = (Info(2).EQ.1) |
---|
28 | C IF Embed3 = .true. THEN the third order embedded formula is used |
---|
29 | C .false. THEN a first order embedded formula is used |
---|
30 | C |
---|
31 | C |
---|
32 | C OUTPUT ARGUMENTS: |
---|
33 | C y = the values of concentrations at TEND. |
---|
34 | C T = equals TEND on output. |
---|
35 | C Info(2) = # of FUN CALLs. |
---|
36 | C Info(3) = # of JAC_SP CALLs. |
---|
37 | C Info(4) = # of accepted steps. |
---|
38 | C Info(5) = # of rejected steps. |
---|
39 | |
---|
40 | INTEGER max_no_steps |
---|
41 | PARAMETER (max_no_steps = 200) |
---|
42 | KPP_REAL Trajectory(NVAR,max_no_steps) |
---|
43 | KPP_REAL StepSize(max_no_steps) |
---|
44 | |
---|
45 | KPP_REAL K1(NVAR), K2(NVAR), K3(NVAR) |
---|
46 | KPP_REAL F1(NVAR), JAC(LU_NONZERO) |
---|
47 | KPP_REAL DFDT(NVAR)(NRAD) |
---|
48 | KPP_REAL Fix(NFIX), Rconst(NREACT) |
---|
49 | KPP_REAL Hmin,Hmax,Hstart,ghinv,uround |
---|
50 | KPP_REAL y(NVAR), Ynew(NVAR) |
---|
51 | KPP_REAL AbsTol(NVAR), RelTol(NVAR) |
---|
52 | KPP_REAL T, Tnext, H, Hold, Tplus |
---|
53 | KPP_REAL ERR, factor, facmax |
---|
54 | KPP_REAL Lambda(NVAR), K11(NVAR), JAC1(LU_NONZERO) |
---|
55 | INTEGER n,nfcn,njac,Naccept,Nreject,i,j |
---|
56 | INTEGER Info(5) |
---|
57 | LOGICAL IsReject, Autonomous, Embed3 |
---|
58 | EXTERNAL FUN, JAC_SP |
---|
59 | |
---|
60 | KPP_REAL gamma, m1, m2, alpha, beta1, beta2, delta, w, e |
---|
61 | KPP_REAL ginv |
---|
62 | c Initialization of counters, etc. |
---|
63 | Autonomous = Info(1) .EQ. 1 |
---|
64 | Embed3 = Info(2) .EQ. 1 |
---|
65 | uround = 1.d-15 |
---|
66 | dround = dsqrt(uround) |
---|
67 | H = DMAX1(Hstart,DMAX1(1.d-8, Hmin)) |
---|
68 | Tplus = T |
---|
69 | IsReject = .false. |
---|
70 | Naccept = 0 |
---|
71 | Nreject = 0 |
---|
72 | Nfcn = 0 |
---|
73 | Njac = 0 |
---|
74 | |
---|
75 | C Method Parameters |
---|
76 | gamma = 1.d0 + 1.d0/sqrt(2.d0) |
---|
77 | a21 = - 1.d0/gamma |
---|
78 | m1 = -3.d0/(2.d0*gamma) |
---|
79 | m2 = -1.d0/(2.d0*gamma) |
---|
80 | c31 = -1.0D0/gamma**2*(1.0D0-7.0D0*gamma+9.0D0*gamma**2) |
---|
81 | & /(-1.0D0+2.0D0*gamma) |
---|
82 | c32 = -1.0D0/gamma**2*(1.0D0-6.0D0*gamma+6.0D0*gamma**2) |
---|
83 | & /(-1.0D0+2.0D0*gamma)/2 |
---|
84 | gamma3 = 0.5D0 - 2*gamma |
---|
85 | d1 = ((-9.0D0*gamma+8.0D0*gamma**2+2.0D0)/gamma**2/ |
---|
86 | & (-1.0D0+2*gamma))/6.0D0 |
---|
87 | d2 = ((-1.0D0+3.0D0*gam)/gamma**2/ |
---|
88 | & (-1.0D0+2.0D0*gamma))/6.0D0 |
---|
89 | d3 = -1.0D0/(3.0D0*gamma) |
---|
90 | |
---|
91 | Trajectory(1:NVAR,1) = Ynew(1) |
---|
92 | |
---|
93 | C === Starting the time loop === |
---|
94 | 10 CONTINUE |
---|
95 | Tplus = T + H |
---|
96 | IF ( Tplus .gt. Tnext ) THEN |
---|
97 | H = Tnext - T |
---|
98 | Tplus = Tnext |
---|
99 | END IF |
---|
100 | |
---|
101 | CALL Jac_SP( Y, Fix, Rconst, JAC ) |
---|
102 | |
---|
103 | Njac = Njac+1 |
---|
104 | ghinv = -1.0d0/(gamma*H) |
---|
105 | DO 20 j=1,NVAR |
---|
106 | JAC(LU_DIAG_V(j)) = JAC(LU_DIAG_V(j)) + ghinv |
---|
107 | 20 CONTINUE |
---|
108 | CALL KppDecomp (NVAR, JAC, ier) |
---|
109 | |
---|
110 | IF (ier.ne.0) THEN |
---|
111 | IF ( H.gt.Hmin) THEN |
---|
112 | H = 5.0d-1*H |
---|
113 | GO TO 10 |
---|
114 | else |
---|
115 | PRINT *,'IER <> 0, H=',H |
---|
116 | STOP |
---|
117 | END IF |
---|
118 | END IF |
---|
119 | |
---|
120 | CALL Fun( Y, Fix, Rconst, F1 ) |
---|
121 | |
---|
122 | C ====== NONAUTONOMOUS CASE =============== |
---|
123 | IF (.not. Autonomous) THEN |
---|
124 | tau = dsign(dround*dmax1( 1.0d-6, dabs(T) ), T) |
---|
125 | CALL Fun( Y, Fix, Rconst, K2 ) |
---|
126 | nfcn=nfcn+1 |
---|
127 | DO 30 j = 1,NVAR |
---|
128 | DFDT(j) = ( K2(j)-F1(j) )/tau |
---|
129 | 30 CONTINUE |
---|
130 | END IF ! .NOT.Autonomous |
---|
131 | |
---|
132 | C ----- STAGE 1 ----- |
---|
133 | DO 40 j = 1,NVAR |
---|
134 | K1(j) = F1(j) |
---|
135 | 40 CONTINUE |
---|
136 | IF (.NOT.Autonomous) THEN |
---|
137 | delta = gamma*H |
---|
138 | DO 45 j = 1,NVAR |
---|
139 | K1(j) = K1(j) + delta*DFDT(j) |
---|
140 | 45 CONTINUE |
---|
141 | END IF ! .NOT.Autonomous |
---|
142 | CALL KppSolve (JAC, K1) |
---|
143 | |
---|
144 | C ----- STAGE 2 ----- |
---|
145 | DO 50 j = 1,NVAR |
---|
146 | Ynew(j) = y(j) + a21*K1(j) |
---|
147 | 50 CONTINUE |
---|
148 | CALL Fun( Ynew, Fix, Rconst, F1 ) |
---|
149 | nfcn=nfcn+1 |
---|
150 | beta = 2.d0/(gamma*H) |
---|
151 | delta = -gamma*H |
---|
152 | DO 55 j = 1,NVAR |
---|
153 | K2(j) = F1(j) + beta*K1(j) |
---|
154 | 55 CONTINUE |
---|
155 | IF (.NOT.Autonomous) THEN |
---|
156 | DO 56 j = 1,NVAR |
---|
157 | K2(j) = K2(j) + delta*DFDT(j) |
---|
158 | 56 CONTINUE |
---|
159 | END IF ! .NOT.Autonomous |
---|
160 | CALL KppSolve (JAC, K2) |
---|
161 | |
---|
162 | C ----- STAGE 3 ----- |
---|
163 | IF (Embed3) THEN |
---|
164 | beta1 = -c31/H |
---|
165 | beta2 = -c32/H |
---|
166 | delta = gamma3*H |
---|
167 | DO 57 j = 1,NVAR |
---|
168 | K3(j) = F1(j) + beta1*K1(j) + beta2*K2(j) |
---|
169 | 57 CONTINUE |
---|
170 | IF (.NOT.Autonomous) THEN |
---|
171 | DO 58 j = 1,NVAR |
---|
172 | K3(j) = K3(j) + delta*DFDT(j) |
---|
173 | 58 CONTINUE |
---|
174 | END IF ! .NOT.Autonomous |
---|
175 | CALL KppSolve (JAC, K3) |
---|
176 | END IF ! Embed3 |
---|
177 | |
---|
178 | |
---|
179 | C ---- The Solution --- |
---|
180 | DO 120 j = 1,NVAR |
---|
181 | Ynew(j) = y(j) + m1*K1(j) + m2*K2(j) |
---|
182 | 120 CONTINUE |
---|
183 | |
---|
184 | |
---|
185 | C ====== Error estimation ======== |
---|
186 | |
---|
187 | ERR=0.d0 |
---|
188 | DO 130 i=1,NVAR |
---|
189 | w = AbsTol(i) + RelTol(i)*DMAX1(DABS(y(i)),DABS(Ynew(i))) |
---|
190 | IF ( Embed3 ) THEN |
---|
191 | e = d1*K1(i) + d2*K2(i) + d3*K3(i) |
---|
192 | ELSE |
---|
193 | e = 1.d0/(2.d0*gamma)*(K1(i)+K2(i)) |
---|
194 | END IF ! Embed3 |
---|
195 | ERR = ERR + ( e/w )**2 |
---|
196 | 130 CONTINUE |
---|
197 | ERR = DMAX1( uround, DSQRT( ERR/NVAR ) ) |
---|
198 | |
---|
199 | C ======= Choose the stepsize =============================== |
---|
200 | |
---|
201 | IF ( Embed3 ) THEN |
---|
202 | elo = 3.0D0 ! estimator local order |
---|
203 | ELSE |
---|
204 | elo = 2.0D0 |
---|
205 | END IF |
---|
206 | factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0)) |
---|
207 | Hnew = DMIN1(Hmax,DMAX1(Hmin, H/factor)) |
---|
208 | Hold = H |
---|
209 | |
---|
210 | C ======= Rejected/Accepted Step ============================ |
---|
211 | |
---|
212 | IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN |
---|
213 | IsReject = .true. |
---|
214 | H = DMIN1(H/10,Hnew) |
---|
215 | Nreject = Nreject+1 |
---|
216 | ELSE |
---|
217 | DO 140 i=1,NVAR |
---|
218 | y(i) = Ynew(i) |
---|
219 | 140 CONTINUE |
---|
220 | T = Tplus |
---|
221 | IF (.NOT.IsReject) THEN |
---|
222 | H = Hnew ! Do not increase stepsize IF previous step was rejected |
---|
223 | END IF |
---|
224 | IsReject = .false. |
---|
225 | Naccept = Naccept+1 |
---|
226 | IF (Naccept+1>max_no_steps) THEN |
---|
227 | PRINT*,'Error in Adjoint Ros2: more steps than allowed' |
---|
228 | STOP |
---|
229 | END IF |
---|
230 | Trajectory(1:NVAR,Naccept+1) = Ynew(1:NVAR) |
---|
231 | StepSize(Naccept) = Hold |
---|
232 | ! CALL TRAJISTORE(y,hold) |
---|
233 | END IF |
---|
234 | C ======= END of the time loop =============================== |
---|
235 | IF ( T .lt. Tnext ) GO TO 10 |
---|
236 | |
---|
237 | C ======= Output Information ================================= |
---|
238 | Info(2) = Nfcn |
---|
239 | Info(3) = Njac |
---|
240 | Info(4) = Naccept |
---|
241 | Info(5) = Nreject |
---|
242 | |
---|
243 | ginv = 1.d0/gamma |
---|
244 | C -- The backwards loop for the CONTINUOUS ADJOINT |
---|
245 | DO istep = Naccept,1,-1 |
---|
246 | |
---|
247 | h = StepSize(istep) |
---|
248 | y(1:NVAR) = Trajectory(1:NVAR,istep+1) |
---|
249 | gHinv = -ginv/H |
---|
250 | |
---|
251 | CALL Jac_SP(Y, Fix, Rconst, JAC) |
---|
252 | JAC1(1:LU_NONZERO)=JAC(1:LU_NONZERO) |
---|
253 | DO j=1,NVAR |
---|
254 | JAC(lu_diag_v(j)) = JAC(lu_diag_v(j)) + gHinv |
---|
255 | END DO |
---|
256 | CALL KppDecomp (NVAR,JAC,ier) |
---|
257 | ccc equivalent to function evaluation in forward integration |
---|
258 | ccc is J^T*Lambda in backward integration |
---|
259 | CALL JacTR_SP_Vec ( JAC1, Lambda, F1) |
---|
260 | |
---|
261 | C ----- STAGE 1 (AUTONOMOUS) ----- |
---|
262 | K11(1:NVAR) = F1(1:NVAR) |
---|
263 | CALL KppSolveTR (JAC,K11,K1) |
---|
264 | C ----- STAGE 2 (AUTONOMOUS) ----- |
---|
265 | y(1:NVAR) = Trajectory(1:NVAR,istep) |
---|
266 | CALL Jac_SP(Y, Fix, Rconst, JAC1) |
---|
267 | Ynew(1:NVAR) = Lambda(1:NVAR) - ginv*K1(1:NVAR) |
---|
268 | CALL JacTR_SP_Vec ( JAC1, Ynew, F1) |
---|
269 | beta = -2.d0*ghinv |
---|
270 | K11(1:NVAR) = F1(1:NVAR) + beta*K1(1:NVAR) |
---|
271 | CALL KppSolveTR (JAC,K11,K2) |
---|
272 | c ---- The solution |
---|
273 | Lambda(1:NVAR) = Lambda(1:NVAR)+m1*K1(1:NVAR)+m2*K2(1:NVAR) |
---|
274 | |
---|
275 | END DO ! istep |
---|
276 | |
---|
277 | |
---|
278 | RETURN |
---|
279 | END |
---|
280 | |
---|
281 | |
---|
282 | SUBROUTINE FUNC_CHEM(N, T, Y, P) |
---|
283 | INCLUDE 'KPP_ROOT_params.h' |
---|
284 | INCLUDE 'KPP_ROOT_global.h' |
---|
285 | INTEGER N |
---|
286 | KPP_REAL T, Told |
---|
287 | KPP_REAL Y(NVAR), P(NVAR) |
---|
288 | Told = TIME |
---|
289 | TIME = T |
---|
290 | CALL Update_SUN() |
---|
291 | CALL Update_RCONST() |
---|
292 | CALL Fun( Y, FIX, RCONST, P ) |
---|
293 | TIME = Told |
---|
294 | RETURN |
---|
295 | END |
---|
296 | |
---|
297 | |
---|
298 | SUBROUTINE JAC_CHEM(N, T, Y, J) |
---|
299 | INCLUDE 'KPP_ROOT_params.h' |
---|
300 | INCLUDE 'KPP_ROOT_global.h' |
---|
301 | INTEGER N |
---|
302 | KPP_REAL Told, T |
---|
303 | KPP_REAL Y(NVAR), J(LU_NONZERO) |
---|
304 | Told = TIME |
---|
305 | TIME = T |
---|
306 | CALL Update_SUN() |
---|
307 | CALL Update_RCONST() |
---|
308 | CALL Jac_SP( Y, FIX, RCONST, J ) |
---|
309 | TIME = Told |
---|
310 | RETURN |
---|
311 | END |
---|
312 | |
---|