Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

ros2_ddm.f @ 2696

Last change on this file since 2696 was 2696, checked in by kanani, 6 years ago
Merge of branch palm4u into trunk
File size: 15.0 KB

Line
1	SUBROUTINE INTEGRATE( NSENSIT, Y, TIN, TOUT )
2
3	INCLUDE 'KPP_ROOT_params.h'
4	INCLUDE 'KPP_ROOT_global.h'
5
6	INTEGER NSENSIT
7	C TIN - Start Time
8	KPP_REAL TIN
9	C TOUT - End Time
10	KPP_REAL TOUT
11	C Y - Concentrations and Sensitivities
12	KPP_REAL Y(NVAR*(NSENSIT+1))
13	C --- Note: Y contains: (1:NVAR) concentrations, followed by
14	C --- (1:NVAR) sensitivities w.r.t. first parameter, followed by
15	C --- etc., followed by
16	C --- (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
17
18	INTEGER INFO(5)
19
20	EXTERNAL FUNC_CHEM, JAC_CHEM
21
22	INFO(1) = Autonomous
23
24	CALL ROS2_DDM(NVAR,NSENSIT,TIN,TOUT,STEPMIN,STEPMAX,
25	+ STEPMIN,Y,ATOL,RTOL,
26	+ Info,FUNC_CHEM,JAC_CHEM)
27
28
29	RETURN
30	END
31
32
33
34
35	SUBROUTINE ROS2_DDM(N,NSENSIT,T,Tnext,Hmin,Hmax,Hstart,
36	+ y,AbsTol,RelTol,
37	+ Info,FUNC_CHEM,JAC_CHEM)
38	IMPLICIT NONE
39	INCLUDE 'KPP_ROOT_params.h'
40	INCLUDE 'KPP_ROOT_global.h'
41	INCLUDE 'KPP_ROOT_sparse.h'
42	C
43	C Ros2 with direct-decoupled calculation of sensitivities
44	C
45	C The global variable DDMTYPE distinguishes between:
46	C DDMTYPE = 0 : sensitivities w.r.t. initial values
47	C DDMTYPE = 1 : sensitivities w.r.t. parameters
48	C
49	C INPUT ARGUMENTS:
50	C y = Vector of: (1:NVAR) concentrations, followed by
51	C (1:NVAR) sensitivities w.r.t. first parameter, followed by
52	C etc., followed by
53	C (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
54	C (y contains initial values at input, final values at output)
55	C [T, Tnext] = the integration interval
56	C Hmin, Hmax = lower and upper bounds for the selected step-size.
57	C Note that for Step = Hmin the current computed
58	C solution is unconditionally accepted by the error
59	C control mechanism.
60	C AbsTol, RelTol = (NVAR) dimensional vectors of
61	C componentwise absolute and relative tolerances.
62	C FUNC_CHEM = name of routine of derivatives. KPP syntax.
63	C See the header below.
64	C JAC_CHEM = name of routine that computes the Jacobian, in
65	C sparse format. KPP syntax. See the header below.
66	C Info(1) = 1 for autonomous system
67	C = 0 for nonautonomous system
68	C
69	C OUTPUT ARGUMENTS:
70	C y = the values of concentrations at TEND.
71	C T = equals TEND on output.
72	C Info(2) = # of FUNC_CHEM calls.
73	C Info(3) = # of JAC_CHEM calls.
74	C Info(4) = # of accepted steps.
75	C Info(5) = # of rejected steps.
76	C
77	C Adrian Sandu, December 2001
78
79
80	INTEGER NSENSIT
81	KPP_REAL y(NVAR(NSENSIT+1)), ynew(NVAR(NSENSIT+1))
82	KPP_REAL K1(NVAR*(NSENSIT+1))
83	KPP_REAL K2(NVAR*(NSENSIT+1))
84	KPP_REAL K3(NVAR)
85	KPP_REAL DFDT(NVAR*(NSENSIT+1))
86	KPP_REAL DFDP(NVARNSENSIT+1), DFDPDT(NVARNSENSIT+1)
87	KPP_REAL DJDP(NVAR*NSENSIT+1)
88	KPP_REAL F1(NVAR), F2(NVAR)
89	KPP_REAL JAC(LU_NONZERO), AJAC(LU_NONZERO)
90	KPP_REAL DJDT(LU_NONZERO)
91	KPP_REAL HESS(NHESS)
92	KPP_REAL Hmin,Hmax,Hnew,Hstart,ghinv,uround
93	KPP_REAL AbsTol(NVAR), RelTol(NVAR)
94	KPP_REAL T, Tnext, H, Hold, Tplus, e
95	KPP_REAL ERR, factor, facmax, dround, elo, tau, gam
96
97	INTEGER n,nfcn,njac,Naccept,Nreject,i,j,ier
98	INTEGER Info(5)
99	LOGICAL IsReject,Autonomous,Embed3
100	EXTERNAL FUNC_CHEM, JAC_CHEM, HESS_CHEM
101
102	LOGICAL negative
103	KPP_REAL gamma, m1, m2, alpha, beta, delta, theta, w
104	KPP_REAL gamma3, d1, d2, d3, beta1, beta2
105	KPP_REAL c31, c32, c34
106
107	c Initialization of counters, etc.
108	Autonomous = Info(1) .EQ. 1
109	Embed3 = Info(2) .EQ. 1
110	uround = 1.d-15
111	dround = 1.0d-7 ! DSQRT(uround)
112	H = DMAX1(1.d-8, Hstart)
113	Tplus = T
114	IsReject = .false.
115	Naccept = 0
116	Nreject = 0
117	Nfcn = 0
118	Njac = 0
119	gamma = 1.d0 + 1.d0/DSQRT(2.0d0)
120	c31 = -1.0D0/gamma*2(1.0D0-7.0D0gamma+9.0D0gamma**2)
121	& /(-1.0D0+2.0D0*gamma)
122	c32 = -1.0D0/gamma*2(1.0D0-6.0D0gamma+6.0D0gamma**2)
123	& /(-1.0D0+2.0D0*gamma)/2
124	gamma3 = 0.5D0 - 2*gamma
125	d1 = ((-9.0D0gamma+8.0D0gamma2+2.0D0)/gamma2/
126	& (-1.0D0+2*gamma))/6.0D0
127	d2 = ((-1.0D0+3.0D0gamma)/gamma*2/
128	& (-1.0D0+2.0D0*gamma))/6.0D0
129	d3 = -1.0D0/(3.0D0*gamma)
130	m1 = -3.d0/(2.d0*gamma)
131	m2 = -1.d0/(2.d0*gamma)
132
133
134	C === Starting the time loop ===
135	10 CONTINUE
136	Tplus = T + H
137	IF ( Tplus .gt. Tnext ) THEN
138	H = Tnext - T
139	Tplus = Tnext
140	END IF
141
142	C Initial Function, Jacobian, and Hessian Values
143	CALL FUNC_CHEM(NVAR, T, y, F1)
144	CALL JAC_CHEM(NVAR, T, y, JAC)
145	CALL HESS_CHEM( NVAR, T, y, HESS )
146	IF (DDMTYPE .EQ. 1) THEN
147	CALL DFUNDPAR(NVAR, NSENSIT, T, y, DFDP)
148	END IF
149
150	C Estimate the time derivatives in non-autonomous case
151	IF (.not. Autonomous) THEN
152	tau = DSIGN(dround*DMAX1( 1.0d0, DABS(T) ), T)
153	CALL FUNC_CHEM(NVAR, T+tau, y, K2)
154	nfcn=nfcn+1
155	CALL JAC_CHEM(NVAR, T+tau, y, AJAC)
156	njac=njac+1
157	DO 20 j = 1,NVAR
158	DFDT(j) = ( K2(j)-F1(j) )/tau
159	20 CONTINUE
160	DO 30 j = 1,LU_NONZERO
161	DJDT(j) = ( AJAC(j)-JAC(j) )/tau
162	30 CONTINUE
163	DO 40 i=1,NSENSIT
164	CALL Jac_SP_Vec (DJDT,y(iNVAR+1),DFDT(iNVAR+1))
165	40 CONTINUE
166	END IF ! .not. Autonomous
167
168	Njac = Njac+1
169	ghinv = - 1.0d0/(gamma*H)
170	DO 50 j=1,LU_NONZERO
171	AJAC(j) = JAC(j)
172	50 CONTINUE
173	DO 60 j=1,NVAR
174	AJAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) + ghinv
175	60 CONTINUE
176	CALL KppDecomp (AJAC, ier)
177
178	IF (ier.ne.0) THEN
179	IF ( H.gt.Hmin) THEN
180	H = 5.0d-1*H
181	go to 10
182	ELSE
183	print *,'IER <> 0, H=',H
184	stop
185	END IF
186	END IF
187
188
189
190
191	C ----- STAGE 1 -----
192	delta = gamma*H
193	DO 70 j = 1,NVAR
194	K1(j) = F1(j)
195	70 CONTINUE
196	IF (.NOT. Autonomous) THEN
197	DO 80 j = 1,NVAR
198	K1(j) = K1(j) + delta*DFDT(j)
199	80 CONTINUE
200	END IF
201	CALL KppSolve (AJAC, K1)
202	C --- If derivative w.r.t. parameters
203	IF (DDMTYPE .EQ. 1) THEN
204	CALL DJACDPAR(NVAR, NSENSIT, T, y, K1(1), DJDP)
205	END IF
206	C --- End of derivative w.r.t. parameters
207	DO 120 i=1,NSENSIT
208	CALL Jac_SP_Vec (JAC,y(iNVAR+1),K1(iNVAR+1))
209	CALL Hess_Vec ( HESS, K1(1), y(i*NVAR+1), F2 )
210	DO 90 j=1,NVAR
211	K1(iNVAR+j) = K1(iNVAR+j) + gHinv*F2(j)
212	90 CONTINUE
213	IF (.NOT. Autonomous) THEN
214	DO 100 j = 1,NVAR
215	K1(iNVAR+j) = K1(iNVAR+j) + deltaDFDT(iNVAR+j)
216	100 CONTINUE
217	END IF
218	C --- If derivative w.r.t. parameters
219	IF (DDMTYPE .EQ. 1) THEN
220	DO 110 j = 1,NVAR
221	K1(iNVAR+j) = K1(iNVAR+j) + DFDP((i-1)*NVAR+j)
222	& + DJDP((i-1)*NVAR+j)
223	110 CONTINUE
224	END IF
225	C --- End of derivative w.r.t. parameters
226	CALL KppSolve (AJAC, K1(i*NVAR+1))
227	120 CONTINUE
228
229	C ----- STAGE 2 -----
230	alpha = - 1.d0/gamma
231	DO 130 j = 1,NVAR*(NSENSIT+1)
232	ynew(j) = y(j) + alpha*K1(j)
233	130 CONTINUE
234	CALL FUNC_CHEM(NVAR, T+H, ynew, F1)
235	IF (DDMTYPE.EQ.1) THEN
236	CALL DFUNDPAR(NVAR, NSENSIT, T+H, ynew, DFDP)
237	END IF
238	nfcn=nfcn+1
239	beta1 = 2.d0/(gamma*H)
240	delta = -gamma*H
241	DO 140 j = 1,NVAR
242	K2(j) = F1(j) + beta1*K1(j)
243	140 CONTINUE
244	IF (.NOT. Autonomous) THEN
245	DO 150 j = 1,NVAR
246	K2(j) = K2(j) + delta*DFDT(j)
247	150 CONTINUE
248	END IF
249	CALL KppSolve (AJAC, K2)
250	C --- If derivative w.r.t. parameters
251	IF (DDMTYPE .EQ. 1) THEN
252	CALL DJACDPAR(NVAR, NSENSIT, T, y, K2(1), DJDP)
253	END IF
254	C --- End of derivative w.r.t. parameters
255
256	CALL JAC_CHEM(NVAR, T+H, Ynew, JAC)
257	njac=njac+1
258	DO 190 i=1,NSENSIT
259	CALL Jac_SP_Vec (JAC,ynew(iNVAR+1),K2(iNVAR+1))
260	CALL Jac_SP_Vec (DJDT,y(i*NVAR+1),F1)
261	CALL Hess_Vec ( HESS, K2(1), y(i*NVAR+1), F2 )
262	DO 160 j = 1,NVAR
263	K2(iNVAR+j) = K2(iNVAR+j) + beta1K1(iNVAR+j)
264	& + gHinv*F2(j)
265	160 CONTINUE
266	IF (.NOT. Autonomous) THEN
267	DO 170 j = 1,NVAR
268	K2(iNVAR+j) = K2(iNVAR+j) + deltaDFDT(iNVAR+j)
269	170 CONTINUE
270	END IF
271	C --- If derivative w.r.t. parameters
272	IF (DDMTYPE .EQ. 1) THEN
273	DO 180 j = 1,NVAR
274	K2(iNVAR+j) = K2(iNVAR+j) + DFDP((i-1)*NVAR+j)
275	& + DJDP((i-1)*NVAR+j)
276	180 CONTINUE
277	END IF
278	C --- End of derivative w.r.t. parameters
279	CALL KppSolve (AJAC, K2(i*NVAR+1))
280	190 CONTINUE
281
282	C ----- STAGE 3 for error control only -----
283	IF (Embed3) THEN
284	beta1 = -c31/H
285	beta2 = -c32/H
286	delta = gamma3*H
287	DO 195 j = 1,NVAR
288	K3(j) = F1(j) + beta1K1(j) + beta2K2(j)
289	195 CONTINUE
290	IF (.NOT. Autonomous) THEN
291	DO 196 j = 1,NVAR
292	K3(j) = K3(j) + delta*DFDT(j)
293	196 CONTINUE
294	END IF
295	CALL KppSolve (AJAC, K3)
296	END IF
297
298	C ---- The Solution ---
299	DO 200 j = 1,NVAR*(NSENSIT+1)
300	ynew(j) = y(j) + m1K1(j) + m2K2(j)
301	200 CONTINUE
302
303
304	C ====== Error estimation for concentrations only; this can be easily adapted to
305	C estimate the sensitivity error too ========
306
307	ERR=0.d0
308	DO 210 i=1,NVAR
309	w = AbsTol(i) + RelTol(i)*DMAX1(DABS(y(i)),DABS(ynew(i)))
310	IF (Embed3) THEN
311	e = d1K1(i) + d2K2(i) + d3*K3(i)
312	ELSE
313	e = (1.d0/(2.d0gamma))(K1(i)+K2(i))
314	END IF
315	ERR = ERR + ( e/w )**2
316	210 CONTINUE
317	ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
318
319	C ======= Choose the stepsize ===============================
320
321	IF (Embed3) THEN
322	elo = 3.0D0 ! estimator local order
323	ELSE
324	elo = 2.0D0
325	END IF
326	factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
327	Hnew = DMIN1(Hmax,DMAX1(Hmin, H/factor))
328
329	C ======= Rejected/Accepted Step ============================
330
331	IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
332	IsReject = .true.
333	H = DMIN1(H/10,Hnew)
334	Nreject = Nreject+1
335	ELSE
336	DO 300 i=1,NVAR*(NSENSIT+1)
337	y(i) = ynew(i)
338	300 CONTINUE
339	T = Tplus
340	IF (.NOT.IsReject) THEN
341	H = Hnew ! Do not increase stepsize if previous step was rejected
342	END IF
343	IsReject = .false.
344	Naccept = Naccept+1
345	END IF
346
347	C ======= End of the time loop ===============================
348	IF ( T .lt. Tnext ) GO TO 10
349
350
351
352	C ======= Output Information =================================
353	Info(2) = Nfcn
354	Info(3) = Njac
355	Info(4) = Naccept
356	Info(5) = Nreject
357
358	RETURN
359	END
360
361
362
363	SUBROUTINE FUNC_CHEM(N, T, Y, P)
364	INCLUDE 'KPP_ROOT_params.h'
365	INCLUDE 'KPP_ROOT_global.h'
366	KPP_REAL T, Told
367	KPP_REAL Y(NVAR), P(NVAR)
368	Told = TIME
369	TIME = T
370	CALL Update_SUN()
371	CALL Update_RCONST()
372	CALL Fun( Y, FIX, RCONST, P )
373	TIME = Told
374	RETURN
375	END
376
377
378	SUBROUTINE DFUNDPAR(N, NSENSIT, T, Y, P)
379	C --- Computes the partial derivatives of FUNC_CHEM w.r.t. parameters
380	INCLUDE 'KPP_ROOT_params.h'
381	INCLUDE 'KPP_ROOT_global.h'
382	C --- NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
383	INTEGER NCOEFF, JCOEFF(NREACT)
384	COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
385
386	KPP_REAL T, Told
387	KPP_REAL Y(NVAR), P(NVAR*NSENSIT)
388	Told = TIME
389	TIME = T
390	CALL Update_SUN()
391	CALL Update_RCONST()
392	C
393	IF (DDMTYPE .EQ. 0) THEN
394	C --- Note: the values below are for sensitivities w.r.t. initial values;
395	C --- they may have to be changed for other applications
396	DO j=1,NSENSIT
397	DO i=1,NVAR
398	P(i+NVAR*(j-1)) = 0.0D0
399	END DO
400	END DO
401	ELSE
402	C --- Example: the call below is for sensitivities w.r.t. rate coefficients;
403	C --- JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
404	C --- w.r.t. which one differentiates
405	CALL dFun_dRcoeff( Y, FIX, NCOEFF, JCOEFF, P )
406	END IF
407	TIME = Told
408	RETURN
409	END
410
411	SUBROUTINE DJACDPAR(N, NSENSIT, T, Y, U, P)
412	C --- Computes the partial derivatives of JAC w.r.t. parameters times user vector U
413	INCLUDE 'KPP_ROOT_params.h'
414	INCLUDE 'KPP_ROOT_global.h'
415	C --- NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
416	INTEGER NCOEFF, JCOEFF(NREACT)
417	COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
418
419	KPP_REAL T, Told
420	KPP_REAL Y(NVAR), U(NVAR)
421	KPP_REAL P(NVAR*NSENSIT)
422	Told = TIME
423	TIME = T
424	CALL Update_SUN()
425	CALL Update_RCONST()
426	C
427	IF (DDMTYPE .EQ. 0) THEN
428	C --- Note: the values below are for sensitivities w.r.t. initial values;
429	C --- they may have to be changed for other applications
430	DO j=1,NSENSIT
431	DO i=1,NVAR
432	P(i+NVAR*(j-1)) = 0.0D0
433	END DO
434	END DO
435	ELSE
436	C --- Example: the call below is for sensitivities w.r.t. rate coefficients;
437	C --- JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
438	C --- w.r.t. which one differentiates
439	CALL dJac_dRcoeff( Y, FIX, U, NCOEFF, JCOEFF, P )
440	END IF
441	TIME = Told
442	RETURN
443	END
444
445	SUBROUTINE JAC_CHEM(N, T, Y, J)
446	INCLUDE 'KPP_ROOT_params.h'
447	INCLUDE 'KPP_ROOT_global.h'
448	INTEGER N
449	KPP_REAL Told, T
450	KPP_REAL Y(NVAR), J(LU_NONZERO)
451	Told = TIME
452	TIME = T
453	CALL Update_SUN()
454	CALL Update_RCONST()
455	CALL Jac_SP( Y, FIX, RCONST, J )
456	TIME = Told
457	RETURN
458	END
459
460
461	SUBROUTINE HESS_CHEM(N, T, Y, HESS)
462	INCLUDE 'KPP_ROOT_params.h'
463	INCLUDE 'KPP_ROOT_global.h'
464	INTEGER N
465	KPP_REAL Told, T
466	KPP_REAL Y(NVAR), HESS(NHESS)
467	Told = TIME
468	TIME = T
469	CALL Update_SUN()
470	CALL Update_RCONST()
471	CALL Hessian( Y, FIX, RCONST, HESS )
472	TIME = Told
473	RETURN
474	END
475
476
477
478
479
480

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int.modified_WCOPY/oldies/ros2_ddm.f @ 2696

Download in other formats: