Home

Context Navigation

Rosenbrock_vec.f90 @ 3561

Last change on this file since 3561 was 3458, checked in by kanani, 6 years ago
Reintegrated fixes/changes from branch chemistry
File size: 42.3 KB

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

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/Rosenbrock_vec.f90 @ 3561

Download in other formats: