Home

Context Navigation

rosenbrock_vec.f90.save @ 2834

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

Line
1	SUBROUTINE Rosenbrock(Y,Tstart,Tend, &
2	AbsTol,RelTol, &
3	RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
4	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5	!
6	! Solves the system y'=F(t,y) using a Rosenbrock method defined by:
7	!
8	! G = 1/(H*gamma(1)) - Jac(t0,Y0)
9	! T_i = t0 + Alpha(i)*H
10	! Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
11	! G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
12	! gamma(i)*dF/dT(t0, Y0)
13	! Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
14	!
15	! For details on Rosenbrock methods and their implementation consult:
16	! E. Hairer and G. Wanner
17	! "Solving ODEs II. Stiff and differential-algebraic problems".
18	! Springer series in computational mathematics, Springer-Verlag, 1996.
19	! The codes contained in the book inspired this implementation.
20	!
21	! (C) Adrian Sandu, August 2004
22	! Virginia Polytechnic Institute and State University
23	! Contact: sandu@cs.vt.edu
24	! This implementation is part of KPP - the Kinetic PreProcessor
25	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
26	!
27	!~~~> INPUT ARGUMENTS:
28	!
29	!- Y(NVAR) = vector of initial conditions (at T=Tstart)
30	!- [Tstart,Tend] = time range of integration
31	! (if Tstart>Tend the integration is performed backwards in time)
32	!- RelTol, AbsTol = user precribed accuracy
33	!- SUBROUTINE Fun( T, Y, Ydot ) = ODE function,
34	! returns Ydot = Y' = F(T,Y)
35	!- SUBROUTINE Jac( T, Y, Jcb ) = Jacobian of the ODE function,
36	! returns Jcb = dFun/dY
37	!- ICNTRL(1:20) = integer inputs parameters
38	!- RCNTRL(1:20) = real inputs parameters
39	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
40	!
41	!~~~> OUTPUT ARGUMENTS:
42	!
43	!- Y(NVAR) -> vector of final states (at T->Tend)
44	!- ISTATUS(1:20) -> integer output parameters
45	!- RSTATUS(1:20) -> real output parameters
46	!- IERR -> job status upon return
47	! success (positive value) or
48	! failure (negative value)
49	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
50	!
51	!~~~> INPUT PARAMETERS:
52	!
53	! Note: For input parameters equal to zero the default values of the
54	! corresponding variables are used.
55	!
56	! ICNTRL(1) = 1: F = F(y) Independent of T (AUTONOMOUS)
57	! = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
58	!
59	! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
60	! = 1: AbsTol, RelTol are scalars
61	!
62	! ICNTRL(3) -> selection of a particular Rosenbrock method
63	! = 0 : default method is Rodas3
64	! = 1 : method is Ros2
65	! = 2 : method is Ros3
66	! = 3 : method is Ros4
67	! = 4 : method is Rodas3
68	! = 5: method is Rodas4
69	!
70	! ICNTRL(4) -> maximum number of integration steps
71	! For ICNTRL(4)=0) the default value of 100000 is used
72	!
73	! RCNTRL(1) -> Hmin, lower bound for the integration step size
74	! It is strongly recommended to keep Hmin = ZERO
75	! RCNTRL(2) -> Hmax, upper bound for the integration step size
76	! RCNTRL(3) -> Hstart, starting value for the integration step size
77	!
78	! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2)
79	! RCNTRL(5) -> FacMax, upper bound on step increase factor (default=6)
80	! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections
81	! (default=0.1)
82	! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller
83	! than the predicted value (default=0.9)
84	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
85	!
86	!~~~> OUTPUT PARAMETERS:
87	!
88	! Note: each call to Rosenbrock adds the current no. of fcn calls
89	! to previous value of ISTATUS(1), and similar for the other params.
90	! Set ISTATUS(1:20) = 0 before call to avoid this accumulation.
91	!
92	! ISTATUS(1) = No. of function calls
93	! ISTATUS(2) = No. of jacobian calls
94	! ISTATUS(3) = No. of steps
95	! ISTATUS(4) = No. of accepted steps
96	! ISTATUS(5) = No. of rejected steps (except at the beginning)
97	! ISTATUS(6) = No. of LU decompositions
98	! ISTATUS(7) = No. of forward/backward substitutions
99	! ISTATUS(8) = No. of singular matrix decompositions
100	!
101	! RSTATUS(1) -> Texit, the time corresponding to the
102	! computed Y upon return
103	! RSTATUS(2) -> Hexit, last accepted step before exit
104	! For multiple restarts, use Hexit as Hstart in the following run
105	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
106
107	IMPLICIT NONE
108
109	!~~~> Arguments
110	REAL(kind=dp), INTENT(INOUT) :: Y(:,:)
111	REAL(kind=dp), INTENT(IN) :: Tstart,Tend
112	REAL(kind=dp), INTENT(IN) :: AbsTol(NVAR),RelTol(NVAR)
113	INTEGER, INTENT(IN) :: ICNTRL(20)
114	REAL(kind=dp), INTENT(IN) :: RCNTRL(20)
115	INTEGER, INTENT(INOUT) :: ISTATUS(20)
116	REAL(kind=dp), INTENT(INOUT) :: RSTATUS(20)
117	INTEGER, INTENT(OUT) :: IERR
118	!~~~> The method parameters
119	INTEGER, PARAMETER :: Smax = 6
120	INTEGER :: Method, ros_S
121	REAL(kind=dp), DIMENSION(Smax) :: ros_M, ros_E, ros_Alpha, ros_Gamma
122	REAL(kind=dp), DIMENSION(Smax*(Smax-1)/2) :: ros_A, ros_C
123	REAL(kind=dp) :: ros_ELO
124	LOGICAL, DIMENSION(Smax) :: ros_NewF
125	CHARACTER(LEN=12) :: ros_Name
126	!~~~> Local variables
127	REAL(kind=dp) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
128	REAL(kind=dp) :: Hmin, Hmax, Hstart, Hexit
129	REAL(kind=dp) :: Texit
130	INTEGER :: i, UplimTol, Max_no_steps
131	LOGICAL :: Autonomous, VectorTol
132	!~~~> Parameters
133	REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp
134	REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
135
136	!~~~> Initialize statistics
137	Nfun = ISTATUS(ifun)
138	Njac = ISTATUS(ijac)
139	Nstp = ISTATUS(istp)
140	Nacc = ISTATUS(iacc)
141	Nrej = ISTATUS(irej)
142	Ndec = ISTATUS(idec)
143	Nsol = ISTATUS(isol)
144	Nsng = ISTATUS(isng)
145
146	!~~~> Autonomous or time dependent ODE. Default is time dependent.
147	Autonomous = .NOT.(ICNTRL(1) == 0)
148
149	!~~~> For Scalar tolerances (ICNTRL(2).NE.0) the code uses AbsTol(1) and RelTol(1)
150	! For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
151	IF (ICNTRL(2) == 0) THEN
152	VectorTol = .TRUE.
153	UplimTol = NVAR
154	ELSE
155	VectorTol = .FALSE.
156	UplimTol = 1
157	END IF
158
159	!~~~> The particular Rosenbrock method chosen
160	IF (ICNTRL(3) == 0) THEN
161	Method = 4
162	ELSEIF ( (ICNTRL(3) >= 1).AND.(ICNTRL(3) <= 5) ) THEN
163	Method = ICNTRL(3)
164	ELSE
165	PRINT * , 'User-selected Rosenbrock method: ICNTRL(3)=', Method
166	CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
167	RETURN
168	END IF
169
170	!~~~> The maximum number of steps admitted
171	IF (ICNTRL(4) == 0) THEN
172	Max_no_steps = 100000 ! 200
173	ELSEIF (ICNTRL(4) > 0) THEN
174	Max_no_steps=ICNTRL(4)
175	ELSE
176	PRINT * ,'User-selected max no. of steps: ICNTRL(4)=',ICNTRL(4)
177	CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR)
178	RETURN
179	END IF
180
181	!~~~> Unit roundoff (1+Roundoff>1)
182	Roundoff = WLAMCH('E')
183
184	!~~~> Lower bound on the step size: (positive value)
185	IF (RCNTRL(1) == ZERO) THEN
186	Hmin = ZERO
187	ELSEIF (RCNTRL(1) > ZERO) THEN
188	Hmin = RCNTRL(1)
189	ELSE
190	PRINT * , 'User-selected Hmin: RCNTRL(1)=', RCNTRL(1)
191	CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
192	RETURN
193	END IF
194	!~~~> Upper bound on the step size: (positive value)
195	IF (RCNTRL(2) == ZERO) THEN
196	Hmax = ABS(Tend-Tstart)
197	ELSEIF (RCNTRL(2) > ZERO) THEN
198	Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart))
199	ELSE
200	PRINT * , 'User-selected Hmax: RCNTRL(2)=', RCNTRL(2)
201	CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
202	RETURN
203	END IF
204	!~~~> Starting step size: (positive value)
205	IF (RCNTRL(3) == ZERO) THEN
206	Hstart = MAX(Hmin,DeltaMin)
207	ELSEIF (RCNTRL(3) > ZERO) THEN
208	Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart))
209	ELSE
210	PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
211	CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
212	RETURN
213	END IF
214	!~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax
215	IF (RCNTRL(4) == ZERO) THEN
216	FacMin = 0.2_dp
217	ELSEIF (RCNTRL(4) > ZERO) THEN
218	FacMin = RCNTRL(4)
219	ELSE
220	PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
221	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
222	RETURN
223	END IF
224	IF (RCNTRL(5) == ZERO) THEN
225	FacMax = 6.0_dp
226	ELSEIF (RCNTRL(5) > ZERO) THEN
227	FacMax = RCNTRL(5)
228	ELSE
229	PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
230	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
231	RETURN
232	END IF
233	!~~~> FacRej: Factor to decrease step after 2 succesive rejections
234	IF (RCNTRL(6) == ZERO) THEN
235	FacRej = 0.1_dp
236	ELSEIF (RCNTRL(6) > ZERO) THEN
237	FacRej = RCNTRL(6)
238	ELSE
239	PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
240	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
241	RETURN
242	END IF
243	!~~~> FacSafe: Safety Factor in the computation of new step size
244	IF (RCNTRL(7) == ZERO) THEN
245	FacSafe = 0.9_dp
246	ELSEIF (RCNTRL(7) > ZERO) THEN
247	FacSafe = RCNTRL(7)
248	ELSE
249	PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
250	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
251	RETURN
252	END IF
253	!~~~> Check if tolerances are reasonable
254	DO i=1,UplimTol
255	IF ( (AbsTol(i) <= ZERO) .OR. (RelTol(i) <= 10.0_dp*Roundoff) &
256	.OR. (RelTol(i) >= 1.0_dp) ) THEN
257	PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
258	PRINT * , ' RelTol(',i,') = ',RelTol(i)
259	CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR)
260	RETURN
261	END IF
262	END DO
263
264
265	!~~~> Initialize the particular Rosenbrock method
266	SELECT CASE (Method)
267	CASE (1)
268	CALL Ros2(ros_S, ros_A, ros_C, ros_M, ros_E, &
269	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
270	CASE (2)
271	CALL Ros3(ros_S, ros_A, ros_C, ros_M, ros_E, &
272	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
273	CASE (3)
274	CALL Ros4(ros_S, ros_A, ros_C, ros_M, ros_E, &
275	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
276	CASE (4)
277	CALL Rodas3(ros_S, ros_A, ros_C, ros_M, ros_E, &
278	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
279	CASE (5)
280	CALL Rodas4(ros_S, ros_A, ros_C, ros_M, ros_E, &
281	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
282	CASE DEFAULT
283	PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=', Method
284	CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
285	RETURN
286	END SELECT
287
288	!~~~> CALL Rosenbrock method
289	CALL ros_Integrator(Y,Tstart,Tend,Texit, &
290	AbsTol, RelTol, &
291	! Rosenbrock method coefficients
292	ros_S, ros_M, ros_E, ros_A, ros_C, &
293	ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
294	! Integration parameters
295	Autonomous, VectorTol, Max_no_steps, &
296	Roundoff, Hmin, Hmax, Hstart, Hexit, &
297	FacMin, FacMax, FacRej, FacSafe, &
298	! Error indicator
299	IERR)
300
301
302	!~~~> Collect run statistics
303	ISTATUS(ifun) = Nfun
304	ISTATUS(ijac) = Njac
305	ISTATUS(istp) = Nstp
306	ISTATUS(iacc) = Nacc
307	ISTATUS(irej) = Nrej
308	ISTATUS(idec) = Ndec
309	ISTATUS(isol) = Nsol
310	ISTATUS(isng) = Nsng
311	!~~~> Last T and H
312	RSTATUS(itexit) = Texit
313	RSTATUS(ihexit) = Hexit
314
315	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
316	CONTAINS ! SUBROUTINES internal to Rosenbrock
317	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
318
319	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
320	SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
321	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
322	! Handles all error messages
323	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
324
325	REAL(kind=dp), INTENT(IN) :: T, H
326	INTEGER, INTENT(IN) :: Code
327	INTEGER, INTENT(OUT) :: IERR
328
329	IERR = Code
330	PRINT * , &
331	'Forced exit from Rosenbrock due to the following error:'
332	IF ((Code>=-8).AND.(Code<=-1)) THEN
333	PRINT *, IERR_NAMES(Code)
334	ELSE
335	PRINT *, 'Unknown Error code: ', Code
336	ENDIF
337
338	PRINT *, "T=", T, "and H=", H
339
340	END SUBROUTINE ros_ErrorMsg
341
342	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
343	SUBROUTINE ros_Integrator (Y, Tstart, Tend, Tout, &
344	AbsTol, RelTol, &
345	!~~~> Rosenbrock method coefficients
346	ros_S, ros_M, ros_E, ros_A, ros_C, &
347	ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
348	!~~~> Integration parameters
349	Autonomous, VectorTol, Max_no_steps, &
350	Roundoff, Hmin, Hmax, Hstart, Hexit, &
351	FacMin, FacMax, FacRej, FacSafe, &
352	!~~~> Error indicator
353	IERR )
354	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
355	! Template for the implementation of a generic Rosenbrock method
356	! defined by ros_S (no of stages)
357	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
358	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
359
360	USE messy_kpp_compress, ONLY: kco_initialize, kco_compress, &
361	kco_finalize, cell_done
362	IMPLICIT NONE
363
364	!~~~> Input: the initial condition at Tstart; Output: the solution at T
365	REAL(kind=dp), INTENT(INOUT) :: Y(:,:)
366	!~~~> Input: integration interval
367	REAL(kind=dp), INTENT(IN) :: Tstart,Tend
368	!~~~> Output: time at which the solution is returned (T=Tend if success)
369	REAL(kind=dp), INTENT(OUT) :: Tout
370	!~~~> Input: tolerances
371	REAL(kind=dp), INTENT(IN) :: AbsTol(NVAR), RelTol(NVAR)
372	!~~~> Input: The Rosenbrock method parameters
373	INTEGER, INTENT(IN) :: ros_S
374	REAL(kind=dp), INTENT(IN) :: ros_M(ros_S), ros_E(ros_S), &
375	ros_Alpha(ros_S), ros_A(ros_S*(ros_S-1)/2), &
376	ros_Gamma(ros_S), ros_C(ros_S*(ros_S-1)/2), ros_ELO
377	LOGICAL, INTENT(IN) :: ros_NewF(ros_S)
378	!~~~> Input: integration parameters
379	LOGICAL, INTENT(IN) :: Autonomous, VectorTol
380	REAL(kind=dp), INTENT(IN) :: Hstart, Hmin, Hmax
381	INTEGER, INTENT(IN) :: Max_no_steps
382	REAL(kind=dp), INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
383	!~~~> Output: last accepted step
384	REAL(kind=dp), INTENT(OUT) :: Hexit
385	!~~~> Output: Error indicator
386	INTEGER, INTENT(OUT) :: IERR
387	! ~~~~ Local variables
388	REAL(kind=dp) :: Ynew(VL_DIM,NVAR), Fcn0(VL_DIM,NVAR), Fcn(VL_DIM,NVAR)
389	REAL(kind=dp) :: K(VL_DIM,NVAR*ros_S)
390	REAL(kind=dp) :: Jac0(VL_DIM,LU_NONZERO), Ghimj(VL_DIM,LU_NONZERO)
391	REAL(kind=dp) :: Yerr(VL_DIM,NVAR)
392
393	!~~~> Local parameters
394	INTEGER :: ioffset, i, j, istage
395	REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp
396	REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
397	!~~~> Locally called functions
398	! REAL(kind=dp) WLAMCH
399	! EXTERNAL WLAMCH
400	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
401	real(kind=dp) :: deltat
402	!kk logical,parameter :: fixed_step=.true.
403
404	integer,save :: icount=0
405
406	! vector compression algorithm
407
408	REAL(kind=dp), dimension(VL_dim) :: T
409	REAL(kind=dp), dimension(VL_dim) :: H, Hnew, HC, HC2, HG, Fac, Tau, Err
410	LOGICAL, dimension(VL_dim) :: RejectLastH, RejectMoreH
411	LOGICAL, dimension(VL_dim) :: Accept_step, error_flag
412	integer :: vl_save
413
414	REAL(kind=dp), dimension(size(FIX,1),size(FIX,2)) :: FIX_lo
415
416	!~~~> Initial preparations
417	T(1:VL) = Tstart
418	Hexit = 0.0_dp
419	H(1:VL) = MIN(Hstart,Hmax)
420	Tout = Tend
421	icount = icount+1
422	where (ABS(H(1:VL)) <= 10.0_dp*Roundoff) H(1:VL) = DeltaMin
423
424	deltat=Tend-Tstart
425
426	RejectLastH(1:VL)=.FALSE.
427	RejectMoreH(1:VL)=.FALSE.
428
429	!~~~> Time loop begins below
430
431	Nstp = 0
432	error_flag(1:VL) = .false.
433	vl_save = VL
434	FIX_lo(1:VL,:) = FIX(1:VL,:)
435	Kacc(1:VL) = 0
436	Krej(1:VL) = 0
437
438	! intialize compression algorithm
439
440	IF (.not. l_fixed_step) then
441	if( return_A ) then
442	call kco_initialize (VL, y, Kacc, Krej, a=a)
443	else
444	call kco_initialize (VL, y, Kacc, Krej)
445	end if
446	end if
447
448	TimeLoop: DO
449
450	IF ( Nstp > Max_no_steps ) THEN ! Too many steps
451	CALL ros_ErrorMsg(-6,T(1),H(1),IERR)
452	! mz_pj_20070514: replace 'call abort()' by
453	! error_flag(1:VL) = .TRUE.
454	! NOTE: Nstp is the maximum number of steps of the remaining current vector
455	call abort() !Vorlaeufig
456	RETURN
457	END IF
458
459	Where ( (T(1:VL)+0.1_dp*H(1:VL)) == T(1:VL)) ! Step size too small
460	error_flag(1:VL) = .TRUE.
461	endwhere
462
463	if( any(error_flag)) then
464	CALL ros_ErrorMsg(-7,T(1),H(1),IERR)
465	call abort() !Vorlaeufig
466	RETURN
467	end if
468
469	!~~~> Limit H if necessary to avoid going beyond Tend
470	!kk Hexit = H
471	H(1:VL) = MIN(H(1:VL),ABS(Tend-T(1:VL)))
472
473	!~~~> Compute the function at current time
474	CALL Fun( Y, FIX_lo, RCONST, Fcn0)
475
476
477	!~~~> Compute the Jacobian at current time
478	CALL Jac_SP( Y, FIX_lo, RCONST, Jac0)
479
480
481	!~~~> Repeat step calculation until current step accepted
482
483	CALL ros_PrepareMatrix(H,ros_Gamma(1), Jac0,Ghimj)
484
485	!~~~> Compute the stages
486	Stage: DO istage = 1, ros_S
487
488	! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR)
489	ioffset = NVAR*(istage-1)
490
491	! For the 1st istage the function has been computed previously
492	IF ( istage == 1 ) THEN
493	Fcn(1:VL,1:NVAR) = Fcn0(1:VL,1:NVAR)
494	! istage>1 and a new function evaluation is needed at the current istage
495	ELSEIF ( ros_NewF(istage) ) THEN
496	select case (istage-1)
497	case (1)
498	Ynew(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
499	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(1- 1)+ 1:NVAR1) &
500	ros_A((istage-1)(istage-2)/2+ 1)
501	case (2)
502	Ynew(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
503	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(2- 1)+ 1:NVAR2) &
504	ros_A((istage-1)(istage-2)/2+ 2)
505	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(2- 1)+ 1:NVAR2) &
506	ros_A((istage-1)(istage-2)/2+ 2)
507	case DEFAULT
508	Ynew(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
509	DO j = 1, istage-1
510	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(j- 1)+ 1:NVARj) &
511	ros_A((istage-1)(istage-2)/2+ j)
512	END DO
513	end select
514
515	CALL Fun( Ynew, FIX_lo, RCONST, Fcn)
516	END IF ! if istage == 1 elseif ros_NewF(istage)
517
518	select case (istage-1)
519	case (1)
520	K(1:VL,ioffset+ 1:ioffset+ NVAR) = Fcn(1:VL,1:NVAR)
521	j = 1
522	HC(1:VL) = ros_C((istage-1)*(istage-2)/2+ j)/(H(1:VL))
523	do i=1,NVAR
524	K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) + K(1:VL,NVAR(j- 1)+ i) HC(1:VL)
525	end do
526	case (2)
527	K(1:VL,ioffset+ 1:ioffset+ NVAR) = Fcn(1:VL,1:NVAR)
528	HC(1:VL) = ros_C((istage-1)*(istage-2)/2+ 1)/(H(1:VL))
529	HC2(1:VL) = ros_C((istage-1)*(istage-2)/2+ 2)/(H(1:VL))
530	do i=1,NVAR
531	K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) + K(1:VL,NVAR(1- 1)+ i) HC(1:VL)
532	K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) + K(1:VL,NVAR(2- 1)+ i) HC2(1:VL)
533	end do
534	case DEFAULT
535	K(1:VL,ioffset+ 1:ioffset+ NVAR) = Fcn(1:VL,1:NVAR)
536	DO j = 1, istage-1
537	HC(1:VL) = ros_C((istage-1)*(istage-2)/2+ j)/(H(1:VL))
538	do i=1,NVAR
539	K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) + K(1:VL,NVAR(j- 1)+ i) HC(1:VL)
540	end do
541	END DO
542	end select
543
544	CALL ros_Solve(Ghimj, K(:,ioffset+1:ioffset+NVAR))
545
546	END DO Stage
547
548
549	if(ros_S == 3) then
550	YNEW(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
551	Yerr(1:VL,1:NVAR) = ZERO
552
553	!~~~> Compute the new solution
554	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(1-1)+1:NVAR1) * ros_m(1)
555	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(2-1)+1:NVAR2) * ros_m(2)
556	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(3-1)+1:NVAR3) * ros_m(3)
557
558	!~~~> Compute the error estimation
559
560	Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + K(1:VL,NVAR(1-1)+1:NVAR1) * ros_E(1)
561	Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + K(1:VL,NVAR(2-1)+1:NVAR2) * ros_E(2)
562	Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + K(1:VL,NVAR(3-1)+1:NVAR3) * ros_E(3)
563	else
564	YNEW(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
565	Yerr(1:VL,1:NVAR) = ZERO
566	DO j=1,ros_S
567	!~~~> Compute the new solution
568	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + &
569	K(1:VL,NVAR(j-1)+1:NVARj) * ros_m(j)
570
571	!~~~> Compute the error estimation
572
573	Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + &
574	K(1:VL,NVAR(j-1)+1:NVARj) * ros_E(j)
575	END DO
576	end if
577
578
579	!~~~> Check the error magnitude and adjust step size
580	Nstp = Nstp+1
581
582	if(l_fixed_step) then
583	Err = 0.0
584	H = t_steps(Nstp)*deltat
585	Hnew = t_steps(Nstp)*deltat
586	else
587	Err(1:VL) = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
588	!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
589	where(Err(1:VL) >= 10.0_dp*Roundoff)
590	Fac(1:VL) = MIN(FacMax,MAX(FacMin,FacSafe/Err(1:VL)**(ONE/ros_ELO)))
591	elsewhere
592	Fac(1:VL) = 1000. !Force large time step
593	endwhere
594	Hnew(1:VL) = H(1:VL)*Fac(1:VL)
595	end if
596
597	Accept_step(1:VL) = (Err(1:VL) <= ONE).OR.(H(1:VL) <= Hmin)
598	where (Accept_step(1:VL))
599	Kacc(1:VL) = Kacc(1:VL)+1
600	T(1:VL) = T(1:VL) + H(1:VL)
601	Hnew(1:VL) = MAX(Hmin,MIN(Hnew(1:VL),Hmax))
602	endwhere
603	do i=1,NVAR
604	where (Accept_step(1:VL))
605	Y(1:VL,i) = MAX(Ynew(1:VL,i),ZERO)
606	endwhere
607	end do
608	where (Accept_step(1:VL) .and. RejectLastH(1:VL))
609	Hnew(1:VL) = MIN(Hnew(1:VL),H(1:VL))
610	endwhere
611	where (Accept_step(1:VL))
612	RejectLastH(1:VL) = .FALSE.
613	RejectMoreH(1:VL) = .FALSE.
614	H(1:VL) = Hnew(1:VL)
615	endwhere
616	!kk ELSE !~~~> Reject step
617	where (.NOT. Accept_step(1:VL) .and. RejectMoreH(1:VL))
618	Hnew(1:VL) = H(1:VL)*FacRej
619	endwhere
620	where (.NOT. Accept_step(1:VL))
621	RejectMoreH(1:VL) = RejectLastH(1:VL)
622	RejectLastH(1:VL) = .TRUE.
623	H(1:VL) = Hnew(1:VL)
624	endwhere
625	where ( .NOT. Accept_step(1:VL) .and. Kacc(1:VL) >= 1)
626	Krej(1:VL) = Krej(1:VL)+1
627	endwhere
628
629
630	if(l_fixed_step) then
631	if(Nstp >= Nfsteps) EXIT
632	else
633	cell_done(1:VL) = ( .NOT. (T(1:VL)-Tend)+Roundoff <= ZERO )
634	call kco_compress (VL, T, H, Hnew, error_flag, y, RCONST, FIX_lo, &
635	RejectLastH, RejectMoreH, Kacc, Krej, a=a)
636	if(VL == 0) EXIT
637	end if
638
639	!kk write(0,*) ' aaa ',icount,Nstp,count(cell_done),count(Accept_step),VL,minval(T)
640	END DO TimeLoop
641
642	!~~~> Succesful exit
643	IERR = 1 !~~~> The integration was successful
644
645	IF (.not. l_fixed_step) then
646	call kco_finalize (y, Kacc, Krej, a=a)
647	end if
648
649	VL = VL_save
650
651	!kk write(0,*) 'bbb ',icount,Nstp,maxval(Kacc),maxval(Krej)
652	!kk if(icount >= 10) call abort()
653
654	END SUBROUTINE ros_Integrator
655
656
657	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
658	FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
659	AbsTol, RelTol, VectorTol )
660	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
661	!~~~> Computes the "scaled norm" of the error vector Yerr
662	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
663	IMPLICIT NONE
664
665	! Input arguments
666	REAL(kind=dp), INTENT(IN) :: Y(:,:), Ynew(:,:), &
667	Yerr(:,:), AbsTol(NVAR), RelTol(NVAR)
668	LOGICAL, INTENT(IN) :: VectorTol
669	! Local variables
670	REAL(kind=dp), dimension(VL) :: Err, Scale, Ymax
671	INTEGER :: i
672	REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp
673	REAL(kind=dp), dimension(VL) :: ros_ErrorNorm
674
675	Err = ZERO
676	DO i=1,NVAR
677	Ymax = MAX( ABS(Y(1:VL,i)),ABS(Ynew(1:VL,i)))
678	IF (VectorTol) THEN
679	Scale = AbsTol(i)+ RelTol(i)*Ymax
680	ELSE
681	Scale = AbsTol(1)+ RelTol(1)*Ymax
682	END IF
683	Err = Err+ (Yerr(1:VL,i)/Scale)**2
684	END DO
685	Err = SQRT(Err/NVAR)
686
687	ros_ErrorNorm = Err
688
689	END FUNCTION ros_ErrorNorm
690
691
692	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
693
694	SUBROUTINE ros_PrepareMatrix ( H, gam, Jac0, Ghimj )
695
696	! --- --- --- --- --- --- --- --- --- --- --- --- ---
697	! Prepares the LHS matrix for stage calculations
698	! 1. Construct Ghimj = 1/(H*ham) - Jac0
699	! "(Gamma H) Inverse Minus Jacobian"
700	! 2. Repeat LU decomposition of Ghimj until successful.
701	! -half the step size if LU decomposition fails and retry
702	! -exit after 5 consecutive fails
703	! --- --- --- --- --- --- --- --- --- --- --- --- ---
704	IMPLICIT NONE
705
706	!~~~> Input arguments
707	REAL(kind=dp), INTENT(IN) :: Jac0(:,:)
708	REAL(kind=dp), INTENT(IN) :: gam
709	!~~~> Output arguments
710	REAL(kind=dp), INTENT(OUT) :: Ghimj(:,:)
711	!~~~> Inout arguments
712	REAL(kind=dp), INTENT(INOUT) :: H(:) ! step size is decreased when LU fails
713	!~~~> Local variables
714	INTEGER :: i, ising, Nconsecutive
715	REAL(kind=dp), dimension(size(H)) :: ghinv
716	REAL(kind=dp), PARAMETER :: ONE = 1.0_dp, HALF = 0.5_dp
717
718	Nconsecutive = 0
719
720
721	!~~~> Construct Ghimj = 1/(H*gam) - Jac0
722	Ghimj(1:VL,1:LU_NONZERO) = -ONE * JAC0(1:VL,1:LU_NONZERO)
723	ghinv(1:VL) = ONE/(H(1:VL)*gam)
724	DO i=1,NVAR
725	Ghimj(1:VL,LU_DIAG(i)) = Ghimj(1:VL,LU_DIAG(i))+ghinv(1:VL)
726	END DO
727
728	!~~~> Compute LU decomposition
729	!CDIR noiexpand
730	CALL KppDecomp( Ghimj, ising )
731
732	Ndec = Ndec + VL
733
734	END SUBROUTINE ros_PrepareMatrix
735
736
737	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
738	SUBROUTINE ros_Solve( A, b )
739	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
740	! Template for the forward/backward substitution (using pre-computed LU decomposition)
741	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
742	IMPLICIT NONE
743	!~~~> Input variables
744	REAL(kind=dp), INTENT(IN) :: A(:,:)
745	!~~~> InOut variables
746	REAL(kind=dp), INTENT(INOUT) :: b(:,:)
747
748	CALL KppSolve( A, b )
749
750	Nsol = Nsol+VL
751
752	END SUBROUTINE ros_Solve
753
754
755
756	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
757	SUBROUTINE Ros2 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
758	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
759	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
760	! --- AN L-STABLE METHOD, 2 stages, order 2
761	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
762
763	IMPLICIT NONE
764
765	INTEGER, PARAMETER :: S=2
766	INTEGER, INTENT(OUT) :: ros_S
767	REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
768	REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
769	REAL(kind=dp), INTENT(OUT) :: ros_ELO
770	LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
771	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
772	DOUBLE PRECISION g
773
774	g = 1.0_dp + 1.0_dp/SQRT(2.0_dp)
775
776	!~~~> Name of the method
777	ros_Name = 'ROS-2'
778	!~~~> Number of stages
779	ros_S = S
780
781	!~~~> The coefficient matrices A and C are strictly lower triangular.
782	! The lower triangular (subdiagonal) elements are stored in row-wise order:
783	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
784	! The general mapping formula is:
785	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
786	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
787
788	ros_A(1) = (1.0_dp)/g
789	ros_C(1) = (-2.0_dp)/g
790	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
791	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
792	ros_NewF(1) = .TRUE.
793	ros_NewF(2) = .TRUE.
794	!~~~> M_i = Coefficients for new step solution
795	ros_M(1)= (3.0_dp)/(2.0_dp*g)
796	ros_M(2)= (1.0_dp)/(2.0_dp*g)
797	! E_i = Coefficients for error estimator
798	ros_E(1) = 1.0_dp/(2.0_dp*g)
799	ros_E(2) = 1.0_dp/(2.0_dp*g)
800	!~~~> ros_ELO = estimator of local order - the minimum between the
801	! main and the embedded scheme orders plus one
802	ros_ELO = 2.0_dp
803	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
804	ros_Alpha(1) = 0.0_dp
805	ros_Alpha(2) = 1.0_dp
806	!~~~> Gamma_i = \sum_j gamma_{i,j}
807	ros_Gamma(1) = g
808	ros_Gamma(2) =-g
809
810	END SUBROUTINE Ros2
811
812
813	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
814	SUBROUTINE Ros3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
815	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
816	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
817	! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
818	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
819
820	IMPLICIT NONE
821
822	INTEGER, PARAMETER :: S=3
823	INTEGER, INTENT(OUT) :: ros_S
824	REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
825	REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
826	REAL(kind=dp), INTENT(OUT) :: ros_ELO
827	LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
828	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
829
830	!~~~> Name of the method
831	ros_Name = 'ROS-3'
832	!~~~> Number of stages
833	ros_S = S
834
835	!~~~> The coefficient matrices A and C are strictly lower triangular.
836	! The lower triangular (subdiagonal) elements are stored in row-wise order:
837	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
838	! The general mapping formula is:
839	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
840	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
841
842	ros_A(1)= 1.0_dp
843	ros_A(2)= 1.0_dp
844	ros_A(3)= 0.0_dp
845
846	ros_C(1) = -0.10156171083877702091975600115545E+01_dp
847	ros_C(2) = 0.40759956452537699824805835358067E+01_dp
848	ros_C(3) = 0.92076794298330791242156818474003E+01_dp
849	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
850	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
851	ros_NewF(1) = .TRUE.
852	ros_NewF(2) = .TRUE.
853	ros_NewF(3) = .FALSE.
854	!~~~> M_i = Coefficients for new step solution
855	ros_M(1) = 0.1E+01_dp
856	ros_M(2) = 0.61697947043828245592553615689730E+01_dp
857	ros_M(3) = -0.42772256543218573326238373806514E+00_dp
858	! E_i = Coefficients for error estimator
859	ros_E(1) = 0.5E+00_dp
860	ros_E(2) = -0.29079558716805469821718236208017E+01_dp
861	ros_E(3) = 0.22354069897811569627360909276199E+00_dp
862	!~~~> ros_ELO = estimator of local order - the minimum between the
863	! main and the embedded scheme orders plus 1
864	ros_ELO = 3.0_dp
865	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
866	ros_Alpha(1)= 0.0E+00_dp
867	ros_Alpha(2)= 0.43586652150845899941601945119356E+00_dp
868	ros_Alpha(3)= 0.43586652150845899941601945119356E+00_dp
869	!~~~> Gamma_i = \sum_j gamma_{i,j}
870	ros_Gamma(1)= 0.43586652150845899941601945119356E+00_dp
871	ros_Gamma(2)= 0.24291996454816804366592249683314E+00_dp
872	ros_Gamma(3)= 0.21851380027664058511513169485832E+01_dp
873
874	END SUBROUTINE Ros3
875
876	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
877
878
879	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
880	SUBROUTINE Ros4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
881	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
882	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
883	! L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
884	! L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
885	!
886	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
887	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
888	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
889	! SPRINGER-VERLAG (1990)
890	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
891
892	IMPLICIT NONE
893
894	INTEGER, PARAMETER :: S=4
895	INTEGER, INTENT(OUT) :: ros_S
896	REAL(kind=dp), DIMENSION(4), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
897	REAL(kind=dp), DIMENSION(6), INTENT(OUT) :: ros_A, ros_C
898	REAL(kind=dp), INTENT(OUT) :: ros_ELO
899	LOGICAL, DIMENSION(4), INTENT(OUT) :: ros_NewF
900	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
901	DOUBLE PRECISION g
902
903	!~~~> Name of the method
904	ros_Name = 'ROS-4'
905	!~~~> Number of stages
906	ros_S = S
907
908	!~~~> The coefficient matrices A and C are strictly lower triangular.
909	! The lower triangular (subdiagonal) elements are stored in row-wise order:
910	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
911	! The general mapping formula is:
912	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
913	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
914
915	ros_A(1) = 0.2000000000000000E+01_dp
916	ros_A(2) = 0.1867943637803922E+01_dp
917	ros_A(3) = 0.2344449711399156E+00_dp
918	ros_A(4) = ros_A(2)
919	ros_A(5) = ros_A(3)
920	ros_A(6) = 0.0_dp
921
922	ros_C(1) =-0.7137615036412310E+01_dp
923	ros_C(2) = 0.2580708087951457E+01_dp
924	ros_C(3) = 0.6515950076447975E+00_dp
925	ros_C(4) =-0.2137148994382534E+01_dp
926	ros_C(5) =-0.3214669691237626E+00_dp
927	ros_C(6) =-0.6949742501781779E+00_dp
928	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
929	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
930	ros_NewF(1) = .TRUE.
931	ros_NewF(2) = .TRUE.
932	ros_NewF(3) = .TRUE.
933	ros_NewF(4) = .FALSE.
934	!~~~> M_i = Coefficients for new step solution
935	ros_M(1) = 0.2255570073418735E+01_dp
936	ros_M(2) = 0.2870493262186792E+00_dp
937	ros_M(3) = 0.4353179431840180E+00_dp
938	ros_M(4) = 0.1093502252409163E+01_dp
939	!~~~> E_i = Coefficients for error estimator
940	ros_E(1) =-0.2815431932141155E+00_dp
941	ros_E(2) =-0.7276199124938920E-01_dp
942	ros_E(3) =-0.1082196201495311E+00_dp
943	ros_E(4) =-0.1093502252409163E+01_dp
944	!~~~> ros_ELO = estimator of local order - the minimum between the
945	! main and the embedded scheme orders plus 1
946	ros_ELO = 4.0_dp
947	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
948	ros_Alpha(1) = 0.0_dp
949	ros_Alpha(2) = 0.1145640000000000E+01_dp
950	ros_Alpha(3) = 0.6552168638155900E+00_dp
951	ros_Alpha(4) = ros_Alpha(3)
952	!~~~> Gamma_i = \sum_j gamma_{i,j}
953	ros_Gamma(1) = 0.5728200000000000E+00_dp
954	ros_Gamma(2) =-0.1769193891319233E+01_dp
955	ros_Gamma(3) = 0.7592633437920482E+00_dp
956	ros_Gamma(4) =-0.1049021087100450E+00_dp
957
958	END SUBROUTINE Ros4
959
960	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
961	SUBROUTINE Rodas3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
962	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
963	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
964	! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
965	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
966
967	IMPLICIT NONE
968
969	INTEGER, PARAMETER :: S=4
970	INTEGER, INTENT(OUT) :: ros_S
971	REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
972	REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
973	REAL(kind=dp), INTENT(OUT) :: ros_ELO
974	LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
975	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
976	DOUBLE PRECISION g
977
978	!~~~> Name of the method
979	ros_Name = 'RODAS-3'
980	!~~~> Number of stages
981	ros_S = S
982
983	!~~~> The coefficient matrices A and C are strictly lower triangular.
984	! The lower triangular (subdiagonal) elements are stored in row-wise order:
985	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
986	! The general mapping formula is:
987	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
988	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
989
990	ros_A(1) = 0.0E+00_dp
991	ros_A(2) = 2.0E+00_dp
992	ros_A(3) = 0.0E+00_dp
993	ros_A(4) = 2.0E+00_dp
994	ros_A(5) = 0.0E+00_dp
995	ros_A(6) = 1.0E+00_dp
996
997	ros_C(1) = 4.0E+00_dp
998	ros_C(2) = 1.0E+00_dp
999	ros_C(3) =-1.0E+00_dp
1000	ros_C(4) = 1.0E+00_dp
1001	ros_C(5) =-1.0E+00_dp
1002	ros_C(6) =-(8.0E+00_dp/3.0E+00_dp)
1003
1004	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1005	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1006	ros_NewF(1) = .TRUE.
1007	ros_NewF(2) = .FALSE.
1008	ros_NewF(3) = .TRUE.
1009	ros_NewF(4) = .TRUE.
1010	!~~~> M_i = Coefficients for new step solution
1011	ros_M(1) = 2.0E+00_dp
1012	ros_M(2) = 0.0E+00_dp
1013	ros_M(3) = 1.0E+00_dp
1014	ros_M(4) = 1.0E+00_dp
1015	!~~~> E_i = Coefficients for error estimator
1016	ros_E(1) = 0.0E+00_dp
1017	ros_E(2) = 0.0E+00_dp
1018	ros_E(3) = 0.0E+00_dp
1019	ros_E(4) = 1.0E+00_dp
1020	!~~~> ros_ELO = estimator of local order - the minimum between the
1021	! main and the embedded scheme orders plus 1
1022	ros_ELO = 3.0E+00_dp
1023	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1024	ros_Alpha(1) = 0.0E+00_dp
1025	ros_Alpha(2) = 0.0E+00_dp
1026	ros_Alpha(3) = 1.0E+00_dp
1027	ros_Alpha(4) = 1.0E+00_dp
1028	!~~~> Gamma_i = \sum_j gamma_{i,j}
1029	ros_Gamma(1) = 0.5E+00_dp
1030	ros_Gamma(2) = 1.5E+00_dp
1031	ros_Gamma(3) = 0.0E+00_dp
1032	ros_Gamma(4) = 0.0E+00_dp
1033
1034	END SUBROUTINE Rodas3
1035
1036	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1037	SUBROUTINE Rodas4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
1038	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
1039	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1040	! STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
1041	!
1042	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
1043	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1044	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1045	! SPRINGER-VERLAG (1996)
1046	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1047
1048	IMPLICIT NONE
1049
1050	INTEGER, PARAMETER :: S=6
1051	INTEGER, INTENT(OUT) :: ros_S
1052	REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
1053	REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
1054	REAL(kind=dp), INTENT(OUT) :: ros_ELO
1055	LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
1056	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
1057	DOUBLE PRECISION g
1058
1059	!~~~> Name of the method
1060	ros_Name = 'RODAS-4'
1061	!~~~> Number of stages
1062	ros_S = 6
1063
1064	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1065	ros_Alpha(1) = 0.000_dp
1066	ros_Alpha(2) = 0.386_dp
1067	ros_Alpha(3) = 0.210_dp
1068	ros_Alpha(4) = 0.630_dp
1069	ros_Alpha(5) = 1.000_dp
1070	ros_Alpha(6) = 1.000_dp
1071
1072	!~~~> Gamma_i = \sum_j gamma_{i,j}
1073	ros_Gamma(1) = 0.2500000000000000E+00_dp
1074	ros_Gamma(2) =-0.1043000000000000E+00_dp
1075	ros_Gamma(3) = 0.1035000000000000E+00_dp
1076	ros_Gamma(4) =-0.3620000000000023E-01_dp
1077	ros_Gamma(5) = 0.0_dp
1078	ros_Gamma(6) = 0.0_dp
1079
1080	!~~~> The coefficient matrices A and C are strictly lower triangular.
1081	! The lower triangular (subdiagonal) elements are stored in row-wise order:
1082	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1083	! The general mapping formula is: A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1084	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1085
1086	ros_A(1) = 0.1544000000000000E+01_dp
1087	ros_A(2) = 0.9466785280815826E+00_dp
1088	ros_A(3) = 0.2557011698983284E+00_dp
1089	ros_A(4) = 0.3314825187068521E+01_dp
1090	ros_A(5) = 0.2896124015972201E+01_dp
1091	ros_A(6) = 0.9986419139977817E+00_dp
1092	ros_A(7) = 0.1221224509226641E+01_dp
1093	ros_A(8) = 0.6019134481288629E+01_dp
1094	ros_A(9) = 0.1253708332932087E+02_dp
1095	ros_A(10) =-0.6878860361058950E+00_dp
1096	ros_A(11) = ros_A(7)
1097	ros_A(12) = ros_A(8)
1098	ros_A(13) = ros_A(9)
1099	ros_A(14) = ros_A(10)
1100	ros_A(15) = 1.0E+00_dp
1101
1102	ros_C(1) =-0.5668800000000000E+01_dp
1103	ros_C(2) =-0.2430093356833875E+01_dp
1104	ros_C(3) =-0.2063599157091915E+00_dp
1105	ros_C(4) =-0.1073529058151375E+00_dp
1106	ros_C(5) =-0.9594562251023355E+01_dp
1107	ros_C(6) =-0.2047028614809616E+02_dp
1108	ros_C(7) = 0.7496443313967647E+01_dp
1109	ros_C(8) =-0.1024680431464352E+02_dp
1110	ros_C(9) =-0.3399990352819905E+02_dp
1111	ros_C(10) = 0.1170890893206160E+02_dp
1112	ros_C(11) = 0.8083246795921522E+01_dp
1113	ros_C(12) =-0.7981132988064893E+01_dp
1114	ros_C(13) =-0.3152159432874371E+02_dp
1115	ros_C(14) = 0.1631930543123136E+02_dp
1116	ros_C(15) =-0.6058818238834054E+01_dp
1117
1118	!~~~> M_i = Coefficients for new step solution
1119	ros_M(1) = ros_A(7)
1120	ros_M(2) = ros_A(8)
1121	ros_M(3) = ros_A(9)
1122	ros_M(4) = ros_A(10)
1123	ros_M(5) = 1.0E+00_dp
1124	ros_M(6) = 1.0E+00_dp
1125
1126	!~~~> E_i = Coefficients for error estimator
1127	ros_E(1) = 0.0E+00_dp
1128	ros_E(2) = 0.0E+00_dp
1129	ros_E(3) = 0.0E+00_dp
1130	ros_E(4) = 0.0E+00_dp
1131	ros_E(5) = 0.0E+00_dp
1132	ros_E(6) = 1.0E+00_dp
1133
1134	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1135	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1136	ros_NewF(1) = .TRUE.
1137	ros_NewF(2) = .TRUE.
1138	ros_NewF(3) = .TRUE.
1139	ros_NewF(4) = .TRUE.
1140	ros_NewF(5) = .TRUE.
1141	ros_NewF(6) = .TRUE.
1142
1143	!~~~> ros_ELO = estimator of local order - the minimum between the
1144	! main and the embedded scheme orders plus 1
1145	ros_ELO = 4.0_dp
1146
1147	END SUBROUTINE Rodas4
1148
1149	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1150	! End of the set of internal Rosenbrock subroutines
1151	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1152	END SUBROUTINE Rosenbrock

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

Download in other formats:

Original Format