Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

rosenbrock.f90 @ 4636

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

Line
1	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
2	! Rosenbrock - Implementation of several Rosenbrock methods: !
3	! * Ros2 !
4	! * Ros3 !
5	! * Ros4 !
6	! * Rodas3 !
7	! * Rodas4 !
8	! By default the code employs the KPP sparse linear algebra routines !
9	! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) !
10	! !
11	! (C) Adrian Sandu, August 2004 !
12	! Virginia Polytechnic Institute and State University !
13	! Contact: sandu@cs.vt.edu !
14	! Revised by Philipp Miehe and Adrian Sandu, May 2006 ! !
15	! This implementation is part of KPP - the Kinetic PreProcessor !
16	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
17
18	MODULE KPP_ROOT_Integrator
19
20	USE KPP_ROOT_Parameters, ONLY: NVAR, NFIX, NSPEC, LU_NONZERO
21	USE KPP_ROOT_Global
22	IMPLICIT NONE
23	PUBLIC
24	SAVE
25
26	!~~~> Statistics on the work performed by the Rosenbrock method
27	INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4, &
28	Nrej=5, Ndec=6, Nsol=7, Nsng=8, &
29	Ntexit=1, Nhexit=2, Nhnew = 3
30
31	CONTAINS
32
33	SUBROUTINE INTEGRATE( TIN, TOUT, &
34	ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )
35
36	IMPLICIT NONE
37
38	KPP_REAL, INTENT(IN) :: TIN ! Start Time
39	KPP_REAL, INTENT(IN) :: TOUT ! End Time
40	! Optional input parameters and statistics
41	INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20)
42	KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20)
43	INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
44	KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
45	INTEGER, INTENT(OUT), OPTIONAL :: IERR_U
46
47	KPP_REAL :: RCNTRL(20), RSTATUS(20)
48	INTEGER :: ICNTRL(20), ISTATUS(20), IERR
49
50	INTEGER, SAVE :: Ntotal = 0
51
52	ICNTRL(:) = 0
53	RCNTRL(:) = 0.0_dp
54	ISTATUS(:) = 0
55	RSTATUS(:) = 0.0_dp
56
57	!~~~> fine-tune the integrator:
58	ICNTRL(1) = 0 ! 0 - non-autonomous, 1 - autonomous
59	ICNTRL(2) = 0 ! 0 - vector tolerances, 1 - scalars
60
61	! If optional parameters are given, and if they are >0,
62	! then they overwrite default settings.
63	IF (PRESENT(ICNTRL_U)) THEN
64	WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)
65	END IF
66	IF (PRESENT(RCNTRL_U)) THEN
67	WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:)
68	END IF
69
70
71	CALL Rosenbrock(NVAR,VAR,TIN,TOUT, &
72	ATOL,RTOL, &
73	RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
74
75	!~~~> Debug option: show no of steps
76	! Ntotal = Ntotal + ISTATUS(Nstp)
77	! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')',' O3=', VAR(ind_O3)
78
79	STEPMIN = RSTATUS(Nhexit)
80	! if optional parameters are given for output they
81	! are updated with the return information
82	IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:)
83	IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(:)
84	IF (PRESENT(IERR_U)) IERR_U = IERR
85
86	END SUBROUTINE INTEGRATE
87
88	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
89	SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, &
90	AbsTol,RelTol, &
91	RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
92	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
93	!
94	! Solves the system y'=F(t,y) using a Rosenbrock method defined by:
95	!
96	! G = 1/(H*gamma(1)) - Jac(t0,Y0)
97	! T_i = t0 + Alpha(i)*H
98	! Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
99	! G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
100	! gamma(i)*dF/dT(t0, Y0)
101	! Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
102	!
103	! For details on Rosenbrock methods and their implementation consult:
104	! E. Hairer and G. Wanner
105	! "Solving ODEs II. Stiff and differential-algebraic problems".
106	! Springer series in computational mathematics, Springer-Verlag, 1996.
107	! The codes contained in the book inspired this implementation.
108	!
109	! (C) Adrian Sandu, August 2004
110	! Virginia Polytechnic Institute and State University
111	! Contact: sandu@cs.vt.edu
112	! Revised by Philipp Miehe and Adrian Sandu, May 2006
113	! This implementation is part of KPP - the Kinetic PreProcessor
114	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
115	!
116	!~~~> INPUT ARGUMENTS:
117	!
118	!- Y(N) = vector of initial conditions (at T=Tstart)
119	!- [Tstart,Tend] = time range of integration
120	! (if Tstart>Tend the integration is performed backwards in time)
121	!- RelTol, AbsTol = user precribed accuracy
122	!- SUBROUTINE Fun( T, Y, Ydot ) = ODE function,
123	! returns Ydot = Y' = F(T,Y)
124	!- SUBROUTINE Jac( T, Y, Jcb ) = Jacobian of the ODE function,
125	! returns Jcb = dFun/dY
126	!- ICNTRL(1:20) = integer inputs parameters
127	!- RCNTRL(1:20) = real inputs parameters
128	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
129	!
130	!~~~> OUTPUT ARGUMENTS:
131	!
132	!- Y(N) -> vector of final states (at T->Tend)
133	!- ISTATUS(1:20) -> integer output parameters
134	!- RSTATUS(1:20) -> real output parameters
135	!- IERR -> job status upon return
136	! success (positive value) or
137	! failure (negative value)
138	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
139	!
140	!~~~> INPUT PARAMETERS:
141	!
142	! Note: For input parameters equal to zero the default values of the
143	! corresponding variables are used.
144	!
145	! ICNTRL(1) = 1: F = F(y) Independent of T (AUTONOMOUS)
146	! = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
147	!
148	! ICNTRL(2) = 0: AbsTol, RelTol are N-dimensional vectors
149	! = 1: AbsTol, RelTol are scalars
150	!
151	! ICNTRL(3) -> selection of a particular Rosenbrock method
152	! = 0 : Rodas3 (default)
153	! = 1 : Ros2
154	! = 2 : Ros3
155	! = 3 : Ros4
156	! = 4 : Rodas3
157	! = 5 : Rodas4
158	!
159	! ICNTRL(4) -> maximum number of integration steps
160	! For ICNTRL(4)=0) the default value of 100000 is used
161	!
162	! RCNTRL(1) -> Hmin, lower bound for the integration step size
163	! It is strongly recommended to keep Hmin = ZERO
164	! RCNTRL(2) -> Hmax, upper bound for the integration step size
165	! RCNTRL(3) -> Hstart, starting value for the integration step size
166	!
167	! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2)
168	! RCNTRL(5) -> FacMax, upper bound on step increase factor (default=6)
169	! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections
170	! (default=0.1)
171	! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller
172	! than the predicted value (default=0.9)
173	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
174	!
175	!
176	! OUTPUT ARGUMENTS:
177	! -----------------
178	!
179	! T -> T value for which the solution has been computed
180	! (after successful return T=Tend).
181	!
182	! Y(N) -> Numerical solution at T
183	!
184	! IDID -> Reports on successfulness upon return:
185	! = 1 for success
186	! < 0 for error (value equals error code)
187	!
188	! ISTATUS(1) -> No. of function calls
189	! ISTATUS(2) -> No. of jacobian calls
190	! ISTATUS(3) -> No. of steps
191	! ISTATUS(4) -> No. of accepted steps
192	! ISTATUS(5) -> No. of rejected steps (except at very beginning)
193	! ISTATUS(6) -> No. of LU decompositions
194	! ISTATUS(7) -> No. of forward/backward substitutions
195	! ISTATUS(8) -> No. of singular matrix decompositions
196	!
197	! RSTATUS(1) -> Texit, the time corresponding to the
198	! computed Y upon return
199	! RSTATUS(2) -> Hexit, last accepted step before exit
200	! RSTATUS(3) -> Hnew, last predicted step (not yet taken)
201	! For multiple restarts, use Hnew as Hstart
202	! in the subsequent run
203	!
204	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
205
206	USE KPP_ROOT_Parameters
207	USE KPP_ROOT_LinearAlgebra
208	IMPLICIT NONE
209
210	!~~~> Arguments
211	INTEGER, INTENT(IN) :: N
212	KPP_REAL, INTENT(INOUT) :: Y(N)
213	KPP_REAL, INTENT(IN) :: Tstart,Tend
214	KPP_REAL, INTENT(IN) :: AbsTol(N),RelTol(N)
215	INTEGER, INTENT(IN) :: ICNTRL(20)
216	KPP_REAL, INTENT(IN) :: RCNTRL(20)
217	INTEGER, INTENT(INOUT) :: ISTATUS(20)
218	KPP_REAL, INTENT(INOUT) :: RSTATUS(20)
219	INTEGER, INTENT(OUT) :: IERR
220	!~~~> Parameters of the Rosenbrock method, up to 6 stages
221	INTEGER :: ros_S, rosMethod
222	INTEGER, PARAMETER :: RS2=1, RS3=2, RS4=3, RD3=4, RD4=5, RG3=6
223	KPP_REAL :: ros_A(15), ros_C(15), ros_M(6), ros_E(6), &
224	ros_Alpha(6), ros_Gamma(6), ros_ELO
225	LOGICAL :: ros_NewF(6)
226	CHARACTER(LEN=12) :: ros_Name
227	!~~~> Local variables
228	KPP_REAL :: Roundoff, FacMin, FacMax, FacRej, FacSafe
229	KPP_REAL :: Hmin, Hmax, Hstart
230	KPP_REAL :: Texit
231	INTEGER :: i, UplimTol, Max_no_steps
232	LOGICAL :: Autonomous, VectorTol
233	!~~~> Parameters
234	KPP_REAL, PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp
235	KPP_REAL, PARAMETER :: DeltaMin = 1.0E-5_dp
236
237	!~~~> Initialize statistics
238	ISTATUS(1:8) = 0
239	RSTATUS(1:3) = ZERO
240
241	!~~~> Autonomous or time dependent ODE. Default is time dependent.
242	Autonomous = .NOT.(ICNTRL(1) == 0)
243
244	!~~~> For Scalar tolerances (ICNTRL(2).NE.0) the code uses AbsTol(1) and RelTol(1)
245	! For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
246	IF (ICNTRL(2) == 0) THEN
247	VectorTol = .TRUE.
248	UplimTol = N
249	ELSE
250	VectorTol = .FALSE.
251	UplimTol = 1
252	END IF
253
254	!~~~> Initialize the particular Rosenbrock method selected
255	SELECT CASE (ICNTRL(3))
256	CASE (1)
257	CALL Ros2
258	CASE (2)
259	CALL Ros3
260	CASE (3)
261	CALL Ros4
262	CASE (0,4)
263	CALL Rodas3
264	CASE (5)
265	CALL Rodas4
266	CASE (6)
267	CALL Rang3
268	CASE DEFAULT
269	PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=',ICNTRL(3)
270	CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
271	RETURN
272	END SELECT
273
274	!~~~> The maximum number of steps admitted
275	IF (ICNTRL(4) == 0) THEN
276	Max_no_steps = 200000
277	ELSEIF (ICNTRL(4) > 0) THEN
278	Max_no_steps=ICNTRL(4)
279	ELSE
280	PRINT * ,'User-selected max no. of steps: ICNTRL(4)=',ICNTRL(4)
281	CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR)
282	RETURN
283	END IF
284
285	!~~~> Unit roundoff (1+Roundoff>1)
286	Roundoff = WLAMCH('E')
287
288	!~~~> Lower bound on the step size: (positive value)
289	IF (RCNTRL(1) == ZERO) THEN
290	Hmin = ZERO
291	ELSEIF (RCNTRL(1) > ZERO) THEN
292	Hmin = RCNTRL(1)
293	ELSE
294	PRINT * , 'User-selected Hmin: RCNTRL(1)=', RCNTRL(1)
295	CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
296	RETURN
297	END IF
298	!~~~> Upper bound on the step size: (positive value)
299	IF (RCNTRL(2) == ZERO) THEN
300	Hmax = ABS(Tend-Tstart)
301	ELSEIF (RCNTRL(2) > ZERO) THEN
302	Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart))
303	ELSE
304	PRINT * , 'User-selected Hmax: RCNTRL(2)=', RCNTRL(2)
305	CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
306	RETURN
307	END IF
308	!~~~> Starting step size: (positive value)
309	IF (RCNTRL(3) == ZERO) THEN
310	Hstart = MAX(Hmin,DeltaMin)
311	ELSEIF (RCNTRL(3) > ZERO) THEN
312	Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart))
313	ELSE
314	PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
315	CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
316	RETURN
317	END IF
318	!~~~> Step size can be changed s.t. FacMin < Hnew/Hold < FacMax
319	IF (RCNTRL(4) == ZERO) THEN
320	FacMin = 0.2_dp
321	ELSEIF (RCNTRL(4) > ZERO) THEN
322	FacMin = RCNTRL(4)
323	ELSE
324	PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
325	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
326	RETURN
327	END IF
328	IF (RCNTRL(5) == ZERO) THEN
329	FacMax = 6.0_dp
330	ELSEIF (RCNTRL(5) > ZERO) THEN
331	FacMax = RCNTRL(5)
332	ELSE
333	PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
334	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
335	RETURN
336	END IF
337	!~~~> FacRej: Factor to decrease step after 2 succesive rejections
338	IF (RCNTRL(6) == ZERO) THEN
339	FacRej = 0.1_dp
340	ELSEIF (RCNTRL(6) > ZERO) THEN
341	FacRej = RCNTRL(6)
342	ELSE
343	PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
344	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
345	RETURN
346	END IF
347	!~~~> FacSafe: Safety Factor in the computation of new step size
348	IF (RCNTRL(7) == ZERO) THEN
349	FacSafe = 0.9_dp
350	ELSEIF (RCNTRL(7) > ZERO) THEN
351	FacSafe = RCNTRL(7)
352	ELSE
353	PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
354	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
355	RETURN
356	END IF
357	!~~~> Check if tolerances are reasonable
358	DO i=1,UplimTol
359	IF ( (AbsTol(i) <= ZERO) .OR. (RelTol(i) <= 10.0_dp*Roundoff) &
360	.OR. (RelTol(i) >= 1.0_dp) ) THEN
361	PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
362	PRINT * , ' RelTol(',i,') = ',RelTol(i)
363	CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR)
364	RETURN
365	END IF
366	END DO
367
368
369	!~~~> CALL Rosenbrock method
370	CALL ros_Integrator(Y, Tstart, Tend, Texit, &
371	AbsTol, RelTol, &
372	! Integration parameters
373	Autonomous, VectorTol, Max_no_steps, &
374	Roundoff, Hmin, Hmax, Hstart, &
375	FacMin, FacMax, FacRej, FacSafe, &
376	! Error indicator
377	IERR)
378
379	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
380	CONTAINS ! SUBROUTINES internal to Rosenbrock
381	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
382
383	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
384	SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
385	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
386	! Handles all error messages
387	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
388
389	KPP_REAL, INTENT(IN) :: T, H
390	INTEGER, INTENT(IN) :: Code
391	INTEGER, INTENT(OUT) :: IERR
392
393	IERR = Code
394	PRINT * , &
395	'Forced exit from Rosenbrock due to the following error:'
396
397	SELECT CASE (Code)
398	CASE (-1)
399	PRINT * , '--> Improper value for maximal no of steps'
400	CASE (-2)
401	PRINT * , '--> Selected Rosenbrock method not implemented'
402	CASE (-3)
403	PRINT * , '--> Hmin/Hmax/Hstart must be positive'
404	CASE (-4)
405	PRINT * , '--> FacMin/FacMax/FacRej must be positive'
406	CASE (-5)
407	PRINT * , '--> Improper tolerance values'
408	CASE (-6)
409	PRINT * , '--> No of steps exceeds maximum bound'
410	CASE (-7)
411	PRINT * , '--> Step size too small: T + 10*H = T', &
412	' or H < Roundoff'
413	CASE (-8)
414	PRINT * , '--> Matrix is repeatedly singular'
415	CASE DEFAULT
416	PRINT *, 'Unknown Error code: ', Code
417	END SELECT
418
419	PRINT *, "T=", T, "and H=", H
420
421	END SUBROUTINE ros_ErrorMsg
422
423	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
424	SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, &
425	AbsTol, RelTol, &
426	!~~~> Integration parameters
427	Autonomous, VectorTol, Max_no_steps, &
428	Roundoff, Hmin, Hmax, Hstart, &
429	FacMin, FacMax, FacRej, FacSafe, &
430	!~~~> Error indicator
431	IERR )
432	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
433	! Template for the implementation of a generic Rosenbrock method
434	! defined by ros_S (no of stages)
435	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
436	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
437
438	IMPLICIT NONE
439
440	!~~~> Input: the initial condition at Tstart; Output: the solution at T
441	KPP_REAL, INTENT(INOUT) :: Y(N)
442	!~~~> Input: integration interval
443	KPP_REAL, INTENT(IN) :: Tstart,Tend
444	!~~~> Output: time at which the solution is returned (T=Tend if success)
445	KPP_REAL, INTENT(OUT) :: T
446	!~~~> Input: tolerances
447	KPP_REAL, INTENT(IN) :: AbsTol(N), RelTol(N)
448	!~~~> Input: integration parameters
449	LOGICAL, INTENT(IN) :: Autonomous, VectorTol
450	KPP_REAL, INTENT(IN) :: Hstart, Hmin, Hmax
451	INTEGER, INTENT(IN) :: Max_no_steps
452	KPP_REAL, INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
453	!~~~> Output: Error indicator
454	INTEGER, INTENT(OUT) :: IERR
455	! ~~~~ Local variables
456	KPP_REAL :: Ynew(N), Fcn0(N), Fcn(N)
457	KPP_REAL :: K(N*ros_S), dFdT(N)
458	#ifdef FULL_ALGEBRA
459	KPP_REAL :: Jac0(N,N), Ghimj(N,N)
460	#else
461	KPP_REAL :: Jac0(LU_NONZERO), Ghimj(LU_NONZERO)
462	#endif
463	KPP_REAL :: H, Hnew, HC, HG, Fac, Tau
464	KPP_REAL :: Err, Yerr(N)
465	INTEGER :: Pivot(N), Direction, ioffset, j, istage
466	LOGICAL :: RejectLastH, RejectMoreH, Singular
467	!~~~> Local parameters
468	KPP_REAL, PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp
469	KPP_REAL, PARAMETER :: DeltaMin = 1.0E-5_dp
470	!~~~> Locally called functions
471	! KPP_REAL WLAMCH
472	! EXTERNAL WLAMCH
473	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
474
475
476	!~~~> Initial preparations
477	T = Tstart
478	RSTATUS(Nhexit) = ZERO
479	H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) )
480	IF (ABS(H) <= 10.0_dp*Roundoff) H = DeltaMin
481
482	IF (Tend >= Tstart) THEN
483	Direction = +1
484	ELSE
485	Direction = -1
486	END IF
487	H = Direction*H
488
489	RejectLastH=.FALSE.
490	RejectMoreH=.FALSE.
491
492	!~~~> Time loop begins below
493
494	TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff <= ZERO) &
495	.OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) )
496
497	IF ( ISTATUS(Nstp) > Max_no_steps ) THEN ! Too many steps
498	CALL ros_ErrorMsg(-6,T,H,IERR)
499	RETURN
500	END IF
501	IF ( ((T+0.1_dp*H) == T).OR.(H <= Roundoff) ) THEN ! Step size too small
502	CALL ros_ErrorMsg(-7,T,H,IERR)
503	RETURN
504	END IF
505
506	!~~~> Limit H if necessary to avoid going beyond Tend
507	H = MIN(H,ABS(Tend-T))
508
509	!~~~> Compute the function at current time
510	CALL FunTemplate(T,Y,Fcn0)
511	ISTATUS(Nfun) = ISTATUS(Nfun) + 1
512
513	!~~~> Compute the function derivative with respect to T
514	IF (.NOT.Autonomous) THEN
515	CALL ros_FunTimeDerivative ( T, Roundoff, Y, &
516	Fcn0, dFdT )
517	END IF
518
519	!~~~> Compute the Jacobian at current time
520	CALL JacTemplate(T,Y,Jac0)
521	ISTATUS(Njac) = ISTATUS(Njac) + 1
522
523	!~~~> Repeat step calculation until current step accepted
524	UntilAccepted: DO
525
526	CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1), &
527	Jac0,Ghimj,Pivot,Singular)
528	IF (Singular) THEN ! More than 5 consecutive failed decompositions
529	CALL ros_ErrorMsg(-8,T,H,IERR)
530	RETURN
531	END IF
532
533	!~~~> Compute the stages
534	Stage: DO istage = 1, ros_S
535
536	! Current istage offset. Current istage vector is K(ioffset+1:ioffset+N)
537	ioffset = N*(istage-1)
538
539	! For the 1st istage the function has been computed previously
540	IF ( istage == 1 ) THEN
541	CALL WCOPY(N,Fcn0,1,Fcn,1)
542	! istage>1 and a new function evaluation is needed at the current istage
543	ELSEIF ( ros_NewF(istage) ) THEN
544	CALL WCOPY(N,Y,1,Ynew,1)
545	DO j = 1, istage-1
546	CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), &
547	K(N(j-1)+1:Nj),1,Ynew,1)
548	END DO
549	Tau = T + ros_Alpha(istage)DirectionH
550	CALL FunTemplate(Tau,Ynew,Fcn)
551	ISTATUS(Nfun) = ISTATUS(Nfun) + 1
552	END IF ! if istage == 1 elseif ros_NewF(istage)
553	CALL WCOPY(N,Fcn,1,K(ioffset+1:ioffset+N),1)
554	DO j = 1, istage-1
555	HC = ros_C((istage-1)(istage-2)/2+j)/(DirectionH)
556	CALL WAXPY(N,HC,K(N(j-1)+1:Nj),1,K(ioffset+1:ioffset+N),1)
557	END DO
558	IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
559	HG = DirectionHros_Gamma(istage)
560	CALL WAXPY(N,HG,dFdT,1,K(ioffset+1:ioffset+N),1)
561	END IF
562	CALL ros_Solve(Ghimj, Pivot, K(ioffset+1:ioffset+N))
563
564	END DO Stage
565
566
567	!~~~> Compute the new solution
568	CALL WCOPY(N,Y,1,Ynew,1)
569	DO j=1,ros_S
570	CALL WAXPY(N,ros_M(j),K(N(j-1)+1:Nj),1,Ynew,1)
571	END DO
572
573	!~~~> Compute the error estimation
574	CALL WSCAL(N,ZERO,Yerr,1)
575	DO j=1,ros_S
576	CALL WAXPY(N,ros_E(j),K(N(j-1)+1:Nj),1,Yerr,1)
577	END DO
578	Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
579
580	!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
581	Fac = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
582	Hnew = H*Fac
583
584	!~~~> Check the error magnitude and adjust step size
585	ISTATUS(Nstp) = ISTATUS(Nstp) + 1
586	IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step
587	ISTATUS(Nacc) = ISTATUS(Nacc) + 1
588	CALL WCOPY(N,Ynew,1,Y,1)
589	T = T + Direction*H
590	Hnew = MAX(Hmin,MIN(Hnew,Hmax))
591	IF (RejectLastH) THEN ! No step size increase after a rejected step
592	Hnew = MIN(Hnew,H)
593	END IF
594	RSTATUS(Nhexit) = H
595	RSTATUS(Nhnew) = Hnew
596	RSTATUS(Ntexit) = T
597	RejectLastH = .FALSE.
598	RejectMoreH = .FALSE.
599	H = Hnew
600	EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
601	ELSE !~~~> Reject step
602	IF (RejectMoreH) THEN
603	Hnew = H*FacRej
604	END IF
605	RejectMoreH = RejectLastH
606	RejectLastH = .TRUE.
607	H = Hnew
608	IF (ISTATUS(Nacc) >= 1) ISTATUS(Nrej) = ISTATUS(Nrej) + 1
609	END IF ! Err <= 1
610
611	END DO UntilAccepted
612
613	END DO TimeLoop
614
615	!~~~> Succesful exit
616	IERR = 1 !~~~> The integration was successful
617
618	END SUBROUTINE ros_Integrator
619
620
621	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
622	KPP_REAL FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
623	AbsTol, RelTol, VectorTol )
624	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
625	!~~~> Computes the "scaled norm" of the error vector Yerr
626	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
627	IMPLICIT NONE
628
629	! Input arguments
630	KPP_REAL, INTENT(IN) :: Y(N), Ynew(N), &
631	Yerr(N), AbsTol(N), RelTol(N)
632	LOGICAL, INTENT(IN) :: VectorTol
633	! Local variables
634	KPP_REAL :: Err, Scale, Ymax
635	INTEGER :: i
636	KPP_REAL, PARAMETER :: ZERO = 0.0_dp
637
638	Err = ZERO
639	DO i=1,N
640	Ymax = MAX(ABS(Y(i)),ABS(Ynew(i)))
641	IF (VectorTol) THEN
642	Scale = AbsTol(i)+RelTol(i)*Ymax
643	ELSE
644	Scale = AbsTol(1)+RelTol(1)*Ymax
645	END IF
646	Err = Err+(Yerr(i)/Scale)**2
647	END DO
648	Err = SQRT(Err/N)
649
650	ros_ErrorNorm = MAX(Err,1.0d-10)
651
652	END FUNCTION ros_ErrorNorm
653
654
655	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
656	SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, &
657	Fcn0, dFdT )
658	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
659	!~~~> The time partial derivative of the function by finite differences
660	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
661	IMPLICIT NONE
662
663	!~~~> Input arguments
664	KPP_REAL, INTENT(IN) :: T, Roundoff, Y(N), Fcn0(N)
665	!~~~> Output arguments
666	KPP_REAL, INTENT(OUT) :: dFdT(N)
667	!~~~> Local variables
668	KPP_REAL :: Delta
669	KPP_REAL, PARAMETER :: ONE = 1.0_dp, DeltaMin = 1.0E-6_dp
670
671	Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
672	CALL FunTemplate(T+Delta,Y,dFdT)
673	ISTATUS(Nfun) = ISTATUS(Nfun) + 1
674	CALL WAXPY(N,(-ONE),Fcn0,1,dFdT,1)
675	CALL WSCAL(N,(ONE/Delta),dFdT,1)
676
677	END SUBROUTINE ros_FunTimeDerivative
678
679
680	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
681	SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, &
682	Jac0, Ghimj, Pivot, Singular )
683	! --- --- --- --- --- --- --- --- --- --- --- --- ---
684	! Prepares the LHS matrix for stage calculations
685	! 1. Construct Ghimj = 1/(H*ham) - Jac0
686	! "(Gamma H) Inverse Minus Jacobian"
687	! 2. Repeat LU decomposition of Ghimj until successful.
688	! -half the step size if LU decomposition fails and retry
689	! -exit after 5 consecutive fails
690	! --- --- --- --- --- --- --- --- --- --- --- --- ---
691	IMPLICIT NONE
692
693	!~~~> Input arguments
694	#ifdef FULL_ALGEBRA
695	KPP_REAL, INTENT(IN) :: Jac0(N,N)
696	#else
697	KPP_REAL, INTENT(IN) :: Jac0(LU_NONZERO)
698	#endif
699	KPP_REAL, INTENT(IN) :: gam
700	INTEGER, INTENT(IN) :: Direction
701	!~~~> Output arguments
702	#ifdef FULL_ALGEBRA
703	KPP_REAL, INTENT(OUT) :: Ghimj(N,N)
704	#else
705	KPP_REAL, INTENT(OUT) :: Ghimj(LU_NONZERO)
706	#endif
707	LOGICAL, INTENT(OUT) :: Singular
708	INTEGER, INTENT(OUT) :: Pivot(N)
709	!~~~> Inout arguments
710	KPP_REAL, INTENT(INOUT) :: H ! step size is decreased when LU fails
711	!~~~> Local variables
712	INTEGER :: i, ising, Nconsecutive
713	KPP_REAL :: ghinv
714	KPP_REAL, PARAMETER :: ONE = 1.0_dp, HALF = 0.5_dp
715
716	Nconsecutive = 0
717	Singular = .TRUE.
718
719	DO WHILE (Singular)
720
721	!~~~> Construct Ghimj = 1/(H*gam) - Jac0
722	#ifdef FULL_ALGEBRA
723	CALL WCOPY(N*N,Jac0,1,Ghimj,1)
724	CALL WSCAL(N*N,(-ONE),Ghimj,1)
725	ghinv = ONE/(DirectionHgam)
726	DO i=1,N
727	Ghimj(i,i) = Ghimj(i,i)+ghinv
728	END DO
729	#else
730	CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1)
731	CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
732	ghinv = ONE/(DirectionHgam)
733	DO i=1,N
734	Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv
735	END DO
736	#endif
737	!~~~> Compute LU decomposition
738	CALL ros_Decomp( Ghimj, Pivot, ising )
739	IF (ising == 0) THEN
740	!~~~> If successful done
741	Singular = .FALSE.
742	ELSE ! ising .ne. 0
743	!~~~> If unsuccessful half the step size; if 5 consecutive fails then return
744	ISTATUS(Nsng) = ISTATUS(Nsng) + 1
745	Nconsecutive = Nconsecutive+1
746	Singular = .TRUE.
747	PRINT*,'Warning: LU Decomposition returned ising = ',ising
748	IF (Nconsecutive <= 5) THEN ! Less than 5 consecutive failed decompositions
749	H = H*HALF
750	ELSE ! More than 5 consecutive failed decompositions
751	RETURN
752	END IF ! Nconsecutive
753	END IF ! ising
754
755	END DO ! WHILE Singular
756
757	END SUBROUTINE ros_PrepareMatrix
758
759
760	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
761	SUBROUTINE ros_Decomp( A, Pivot, ising )
762	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
763	! Template for the LU decomposition
764	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
765	IMPLICIT NONE
766	!~~~> Inout variables
767	KPP_REAL, INTENT(INOUT) :: A(LU_NONZERO)
768	!~~~> Output variables
769	INTEGER, INTENT(OUT) :: Pivot(N), ising
770
771	#ifdef FULL_ALGEBRA
772	CALL DGETRF( N, N, A, N, Pivot, ising )
773	#else
774	CALL KppDecomp ( A, ising )
775	Pivot(1) = 1
776	#endif
777	ISTATUS(Ndec) = ISTATUS(Ndec) + 1
778
779	END SUBROUTINE ros_Decomp
780
781
782	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
783	SUBROUTINE ros_Solve( A, Pivot, b )
784	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
785	! Template for the forward/backward substitution (using pre-computed LU decomposition)
786	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
787	IMPLICIT NONE
788	!~~~> Input variables
789	KPP_REAL, INTENT(IN) :: A(LU_NONZERO)
790	INTEGER, INTENT(IN) :: Pivot(N)
791	!~~~> InOut variables
792	KPP_REAL, INTENT(INOUT) :: b(N)
793
794	#ifdef FULL_ALGEBRA
795	CALL DGETRS( 'N', N , 1, A, N, Pivot, b, N, 0 )
796	#else
797	CALL KppSolve( A, b )
798	#endif
799
800	ISTATUS(Nsol) = ISTATUS(Nsol) + 1
801
802	END SUBROUTINE ros_Solve
803
804
805
806	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
807	SUBROUTINE Ros2
808	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
809	! --- AN L-STABLE METHOD, 2 stages, order 2
810	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
811
812	IMPLICIT NONE
813	DOUBLE PRECISION g
814
815	g = 1.0_dp + 1.0_dp/SQRT(2.0_dp)
816	rosMethod = RS2
817	!~~~> Name of the method
818	ros_Name = 'ROS-2'
819	!~~~> Number of stages
820	ros_S = 2
821
822	!~~~> The coefficient matrices A and C are strictly lower triangular.
823	! The lower triangular (subdiagonal) elements are stored in row-wise order:
824	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
825	! The general mapping formula is:
826	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
827	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
828
829	ros_A(1) = (1.0_dp)/g
830	ros_C(1) = (-2.0_dp)/g
831	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
832	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
833	ros_NewF(1) = .TRUE.
834	ros_NewF(2) = .TRUE.
835	!~~~> M_i = Coefficients for new step solution
836	ros_M(1)= (3.0_dp)/(2.0_dp*g)
837	ros_M(2)= (1.0_dp)/(2.0_dp*g)
838	! E_i = Coefficients for error estimator
839	ros_E(1) = 1.0_dp/(2.0_dp*g)
840	ros_E(2) = 1.0_dp/(2.0_dp*g)
841	!~~~> ros_ELO = estimator of local order - the minimum between the
842	! main and the embedded scheme orders plus one
843	ros_ELO = 2.0_dp
844	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
845	ros_Alpha(1) = 0.0_dp
846	ros_Alpha(2) = 1.0_dp
847	!~~~> Gamma_i = \sum_j gamma_{i,j}
848	ros_Gamma(1) = g
849	ros_Gamma(2) =-g
850
851	END SUBROUTINE Ros2
852
853
854	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
855	SUBROUTINE Ros3
856	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
857	! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
858	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
859
860	IMPLICIT NONE
861	rosMethod = RS3
862	!~~~> Name of the method
863	ros_Name = 'ROS-3'
864	!~~~> Number of stages
865	ros_S = 3
866
867	!~~~> The coefficient matrices A and C are strictly lower triangular.
868	! The lower triangular (subdiagonal) elements are stored in row-wise order:
869	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
870	! The general mapping formula is:
871	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
872	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
873
874	ros_A(1)= 1.0_dp
875	ros_A(2)= 1.0_dp
876	ros_A(3)= 0.0_dp
877
878	ros_C(1) = -0.10156171083877702091975600115545E+01_dp
879	ros_C(2) = 0.40759956452537699824805835358067E+01_dp
880	ros_C(3) = 0.92076794298330791242156818474003E+01_dp
881	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
882	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
883	ros_NewF(1) = .TRUE.
884	ros_NewF(2) = .TRUE.
885	ros_NewF(3) = .FALSE.
886	!~~~> M_i = Coefficients for new step solution
887	ros_M(1) = 0.1E+01_dp
888	ros_M(2) = 0.61697947043828245592553615689730E+01_dp
889	ros_M(3) = -0.42772256543218573326238373806514_dp
890	! E_i = Coefficients for error estimator
891	ros_E(1) = 0.5_dp
892	ros_E(2) = -0.29079558716805469821718236208017E+01_dp
893	ros_E(3) = 0.22354069897811569627360909276199_dp
894	!~~~> ros_ELO = estimator of local order - the minimum between the
895	! main and the embedded scheme orders plus 1
896	ros_ELO = 3.0_dp
897	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
898	ros_Alpha(1)= 0.0_dp
899	ros_Alpha(2)= 0.43586652150845899941601945119356_dp
900	ros_Alpha(3)= 0.43586652150845899941601945119356_dp
901	!~~~> Gamma_i = \sum_j gamma_{i,j}
902	ros_Gamma(1)= 0.43586652150845899941601945119356_dp
903	ros_Gamma(2)= 0.24291996454816804366592249683314_dp
904	ros_Gamma(3)= 0.21851380027664058511513169485832E+01_dp
905
906	END SUBROUTINE Ros3
907
908	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
909
910
911	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
912	SUBROUTINE Ros4
913	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
914	! L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
915	! L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
916	!
917	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
918	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
919	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
920	! SPRINGER-VERLAG (1990)
921	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
922
923	IMPLICIT NONE
924
925	rosMethod = RS4
926	!~~~> Name of the method
927	ros_Name = 'ROS-4'
928	!~~~> Number of stages
929	ros_S = 4
930
931	!~~~> The coefficient matrices A and C are strictly lower triangular.
932	! The lower triangular (subdiagonal) elements are stored in row-wise order:
933	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
934	! The general mapping formula is:
935	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
936	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
937
938	ros_A(1) = 0.2000000000000000E+01_dp
939	ros_A(2) = 0.1867943637803922E+01_dp
940	ros_A(3) = 0.2344449711399156_dp
941	ros_A(4) = ros_A(2)
942	ros_A(5) = ros_A(3)
943	ros_A(6) = 0.0_dp
944
945	ros_C(1) =-0.7137615036412310E+01_dp
946	ros_C(2) = 0.2580708087951457E+01_dp
947	ros_C(3) = 0.6515950076447975_dp
948	ros_C(4) =-0.2137148994382534E+01_dp
949	ros_C(5) =-0.3214669691237626_dp
950	ros_C(6) =-0.6949742501781779_dp
951	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
952	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
953	ros_NewF(1) = .TRUE.
954	ros_NewF(2) = .TRUE.
955	ros_NewF(3) = .TRUE.
956	ros_NewF(4) = .FALSE.
957	!~~~> M_i = Coefficients for new step solution
958	ros_M(1) = 0.2255570073418735E+01_dp
959	ros_M(2) = 0.2870493262186792_dp
960	ros_M(3) = 0.4353179431840180_dp
961	ros_M(4) = 0.1093502252409163E+01_dp
962	!~~~> E_i = Coefficients for error estimator
963	ros_E(1) =-0.2815431932141155_dp
964	ros_E(2) =-0.7276199124938920E-01_dp
965	ros_E(3) =-0.1082196201495311_dp
966	ros_E(4) =-0.1093502252409163E+01_dp
967	!~~~> ros_ELO = estimator of local order - the minimum between the
968	! main and the embedded scheme orders plus 1
969	ros_ELO = 4.0_dp
970	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
971	ros_Alpha(1) = 0.0_dp
972	ros_Alpha(2) = 0.1145640000000000E+01_dp
973	ros_Alpha(3) = 0.6552168638155900_dp
974	ros_Alpha(4) = ros_Alpha(3)
975	!~~~> Gamma_i = \sum_j gamma_{i,j}
976	ros_Gamma(1) = 0.5728200000000000_dp
977	ros_Gamma(2) =-0.1769193891319233E+01_dp
978	ros_Gamma(3) = 0.7592633437920482_dp
979	ros_Gamma(4) =-0.1049021087100450_dp
980
981	END SUBROUTINE Ros4
982
983	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
984	SUBROUTINE Rodas3
985	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
986	! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
987	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
988
989	IMPLICIT NONE
990
991	rosMethod = RD3
992	!~~~> Name of the method
993	ros_Name = 'RODAS-3'
994	!~~~> Number of stages
995	ros_S = 4
996
997	!~~~> The coefficient matrices A and C are strictly lower triangular.
998	! The lower triangular (subdiagonal) elements are stored in row-wise order:
999	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1000	! The general mapping formula is:
1001	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1002	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1003
1004	ros_A(1) = 0.0_dp
1005	ros_A(2) = 2.0_dp
1006	ros_A(3) = 0.0_dp
1007	ros_A(4) = 2.0_dp
1008	ros_A(5) = 0.0_dp
1009	ros_A(6) = 1.0_dp
1010
1011	ros_C(1) = 4.0_dp
1012	ros_C(2) = 1.0_dp
1013	ros_C(3) =-1.0_dp
1014	ros_C(4) = 1.0_dp
1015	ros_C(5) =-1.0_dp
1016	ros_C(6) =-(8.0_dp/3.0_dp)
1017
1018	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1019	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1020	ros_NewF(1) = .TRUE.
1021	ros_NewF(2) = .FALSE.
1022	ros_NewF(3) = .TRUE.
1023	ros_NewF(4) = .TRUE.
1024	!~~~> M_i = Coefficients for new step solution
1025	ros_M(1) = 2.0_dp
1026	ros_M(2) = 0.0_dp
1027	ros_M(3) = 1.0_dp
1028	ros_M(4) = 1.0_dp
1029	!~~~> E_i = Coefficients for error estimator
1030	ros_E(1) = 0.0_dp
1031	ros_E(2) = 0.0_dp
1032	ros_E(3) = 0.0_dp
1033	ros_E(4) = 1.0_dp
1034	!~~~> ros_ELO = estimator of local order - the minimum between the
1035	! main and the embedded scheme orders plus 1
1036	ros_ELO = 3.0_dp
1037	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1038	ros_Alpha(1) = 0.0_dp
1039	ros_Alpha(2) = 0.0_dp
1040	ros_Alpha(3) = 1.0_dp
1041	ros_Alpha(4) = 1.0_dp
1042	!~~~> Gamma_i = \sum_j gamma_{i,j}
1043	ros_Gamma(1) = 0.5_dp
1044	ros_Gamma(2) = 1.5_dp
1045	ros_Gamma(3) = 0.0_dp
1046	ros_Gamma(4) = 0.0_dp
1047
1048	END SUBROUTINE Rodas3
1049
1050	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1051	SUBROUTINE Rodas4
1052	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1053	! STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
1054	!
1055	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
1056	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1057	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1058	! SPRINGER-VERLAG (1996)
1059	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1060
1061	IMPLICIT NONE
1062
1063	rosMethod = RD4
1064	!~~~> Name of the method
1065	ros_Name = 'RODAS-4'
1066	!~~~> Number of stages
1067	ros_S = 6
1068
1069	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1070	ros_Alpha(1) = 0.000_dp
1071	ros_Alpha(2) = 0.386_dp
1072	ros_Alpha(3) = 0.210_dp
1073	ros_Alpha(4) = 0.630_dp
1074	ros_Alpha(5) = 1.000_dp
1075	ros_Alpha(6) = 1.000_dp
1076
1077	!~~~> Gamma_i = \sum_j gamma_{i,j}
1078	ros_Gamma(1) = 0.2500000000000000_dp
1079	ros_Gamma(2) =-0.1043000000000000_dp
1080	ros_Gamma(3) = 0.1035000000000000_dp
1081	ros_Gamma(4) =-0.3620000000000023E-01_dp
1082	ros_Gamma(5) = 0.0_dp
1083	ros_Gamma(6) = 0.0_dp
1084
1085	!~~~> The coefficient matrices A and C are strictly lower triangular.
1086	! The lower triangular (subdiagonal) elements are stored in row-wise order:
1087	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1088	! The general mapping formula is: A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1089	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1090
1091	ros_A(1) = 0.1544000000000000E+01_dp
1092	ros_A(2) = 0.9466785280815826_dp
1093	ros_A(3) = 0.2557011698983284_dp
1094	ros_A(4) = 0.3314825187068521E+01_dp
1095	ros_A(5) = 0.2896124015972201E+01_dp
1096	ros_A(6) = 0.9986419139977817_dp
1097	ros_A(7) = 0.1221224509226641E+01_dp
1098	ros_A(8) = 0.6019134481288629E+01_dp
1099	ros_A(9) = 0.1253708332932087E+02_dp
1100	ros_A(10) =-0.6878860361058950_dp
1101	ros_A(11) = ros_A(7)
1102	ros_A(12) = ros_A(8)
1103	ros_A(13) = ros_A(9)
1104	ros_A(14) = ros_A(10)
1105	ros_A(15) = 1.0_dp
1106
1107	ros_C(1) =-0.5668800000000000E+01_dp
1108	ros_C(2) =-0.2430093356833875E+01_dp
1109	ros_C(3) =-0.2063599157091915_dp
1110	ros_C(4) =-0.1073529058151375_dp
1111	ros_C(5) =-0.9594562251023355E+01_dp
1112	ros_C(6) =-0.2047028614809616E+02_dp
1113	ros_C(7) = 0.7496443313967647E+01_dp
1114	ros_C(8) =-0.1024680431464352E+02_dp
1115	ros_C(9) =-0.3399990352819905E+02_dp
1116	ros_C(10) = 0.1170890893206160E+02_dp
1117	ros_C(11) = 0.8083246795921522E+01_dp
1118	ros_C(12) =-0.7981132988064893E+01_dp
1119	ros_C(13) =-0.3152159432874371E+02_dp
1120	ros_C(14) = 0.1631930543123136E+02_dp
1121	ros_C(15) =-0.6058818238834054E+01_dp
1122
1123	!~~~> M_i = Coefficients for new step solution
1124	ros_M(1) = ros_A(7)
1125	ros_M(2) = ros_A(8)
1126	ros_M(3) = ros_A(9)
1127	ros_M(4) = ros_A(10)
1128	ros_M(5) = 1.0_dp
1129	ros_M(6) = 1.0_dp
1130
1131	!~~~> E_i = Coefficients for error estimator
1132	ros_E(1) = 0.0_dp
1133	ros_E(2) = 0.0_dp
1134	ros_E(3) = 0.0_dp
1135	ros_E(4) = 0.0_dp
1136	ros_E(5) = 0.0_dp
1137	ros_E(6) = 1.0_dp
1138
1139	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1140	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1141	ros_NewF(1) = .TRUE.
1142	ros_NewF(2) = .TRUE.
1143	ros_NewF(3) = .TRUE.
1144	ros_NewF(4) = .TRUE.
1145	ros_NewF(5) = .TRUE.
1146	ros_NewF(6) = .TRUE.
1147
1148	!~~~> ros_ELO = estimator of local order - the minimum between the
1149	! main and the embedded scheme orders plus 1
1150	ros_ELO = 4.0_dp
1151
1152	END SUBROUTINE Rodas4
1153	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1154	SUBROUTINE Rang3
1155	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1156	! STIFFLY-STABLE W METHOD OF ORDER 3, WITH 4 STAGES
1157	!
1158	! J. RANG and L. ANGERMANN
1159	! NEW ROSENBROCK W-METHODS OF ORDER 3
1160	! FOR PARTIAL DIFFERENTIAL ALGEBRAIC
1161	! EQUATIONS OF INDEX 1
1162	! BIT Numerical Mathematics (2005) 45: 761-787
1163	! DOI: 10.1007/s10543-005-0035-y
1164	! Table 4.1-4.2
1165	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1166
1167	IMPLICIT NONE
1168
1169	rosMethod = RG3
1170	!~~~> Name of the method
1171	ros_Name = 'RANG-3'
1172	!~~~> Number of stages
1173	ros_S = 4
1174
1175	ros_A(1) = 5.09052051067020d+00;
1176	ros_A(2) = 5.09052051067020d+00;
1177	ros_A(3) = 0.0d0;
1178	ros_A(4) = 4.97628111010787d+00;
1179	ros_A(5) = 2.77268164715849d-02;
1180	ros_A(6) = 2.29428036027904d-01;
1181
1182	ros_C(1) = -1.16790812312283d+01;
1183	ros_C(2) = -1.64057326467367d+01;
1184	ros_C(3) = -2.77268164715850d-01;
1185	ros_C(4) = -8.38103960500476d+00;
1186	ros_C(5) = -8.48328409199343d-01;
1187	ros_C(6) = 2.87009860433106d-01;
1188
1189	ros_M(1) = 5.22582761233094d+00;
1190	ros_M(2) = -5.56971148154165d-01;
1191	ros_M(3) = 3.57979469353645d-01;
1192	ros_M(4) = 1.72337398521064d+00;
1193
1194	ros_E(1) = -5.16845212784040d+00;
1195	ros_E(2) = -1.26351942603842d+00;
1196	ros_E(3) = -1.11022302462516d-16;
1197	ros_E(4) = 2.22044604925031d-16;
1198
1199	ros_Alpha(1) = 0.0d00;
1200	ros_Alpha(2) = 2.21878746765329d+00;
1201	ros_Alpha(3) = 2.21878746765329d+00;
1202	ros_Alpha(4) = 1.55392337535788d+00;
1203
1204	ros_Gamma(1) = 4.35866521508459d-01;
1205	ros_Gamma(2) = -1.78292094614483d+00;
1206	ros_Gamma(3) = -2.46541900496934d+00;
1207	ros_Gamma(4) = -8.05529997906370d-01;
1208
1209
1210	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1211	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1212	ros_NewF(1) = .TRUE.
1213	ros_NewF(2) = .TRUE.
1214	ros_NewF(3) = .TRUE.
1215	ros_NewF(4) = .TRUE.
1216
1217	!~~~> ros_ELO = estimator of local order - the minimum between the
1218	! main and the embedded scheme orders plus 1
1219	ros_ELO = 3.0_dp
1220
1221	END SUBROUTINE Rang3
1222	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1223
1224	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1225	! End of the set of internal Rosenbrock subroutines
1226	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1227	END SUBROUTINE Rosenbrock
1228	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1229
1230
1231	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1232	SUBROUTINE FunTemplate( T, Y, Ydot )
1233	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1234	! Template for the ODE function call.
1235	! Updates the rate coefficients (and possibly the fixed species) at each call
1236	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1237	USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO
1238	USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME
1239	USE KPP_ROOT_Function, ONLY: Fun
1240	USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST
1241	!~~~> Input variables
1242	KPP_REAL :: T, Y(NVAR)
1243	!~~~> Output variables
1244	KPP_REAL :: Ydot(NVAR)
1245	!~~~> Local variables
1246	KPP_REAL :: Told
1247
1248	Told = TIME
1249	TIME = T
1250	CALL Update_SUN()
1251	CALL Update_RCONST()
1252	CALL Fun( Y, FIX, RCONST, Ydot )
1253	TIME = Told
1254
1255	END SUBROUTINE FunTemplate
1256
1257
1258	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1259	SUBROUTINE JacTemplate( T, Y, Jcb )
1260	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1261	! Template for the ODE Jacobian call.
1262	! Updates the rate coefficients (and possibly the fixed species) at each call
1263	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1264	USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO
1265	USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME
1266	USE KPP_ROOT_Jacobian, ONLY: Jac_SP, LU_IROW, LU_ICOL
1267	USE KPP_ROOT_LinearAlgebra
1268	USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST
1269	!~~~> Input variables
1270	KPP_REAL :: T, Y(NVAR)
1271	!~~~> Output variables
1272	#ifdef FULL_ALGEBRA
1273	KPP_REAL :: JV(LU_NONZERO), Jcb(NVAR,NVAR)
1274	#else
1275	KPP_REAL :: Jcb(LU_NONZERO)
1276	#endif
1277	!~~~> Local variables
1278	KPP_REAL :: Told
1279	#ifdef FULL_ALGEBRA
1280	INTEGER :: i, j
1281	#endif
1282
1283	Told = TIME
1284	TIME = T
1285	CALL Update_SUN()
1286	CALL Update_RCONST()
1287	#ifdef FULL_ALGEBRA
1288	CALL Jac_SP(Y, FIX, RCONST, JV)
1289	DO j=1,NVAR
1290	DO i=1,NVAR
1291	Jcb(i,j) = 0.0_dp
1292	END DO
1293	END DO
1294	DO i=1,LU_NONZERO
1295	Jcb(LU_IROW(i),LU_ICOL(i)) = JV(i)
1296	END DO
1297	#else
1298	CALL Jac_SP( Y, FIX, RCONST, Jcb )
1299	#endif
1300	TIME = Told
1301
1302	END SUBROUTINE JacTemplate
1303
1304	END MODULE KPP_ROOT_Integrator
1305
1306
1307
1308

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int.modified_WCOPY/rosenbrock.f90 @ 4636

Download in other formats: