Home

Context Navigation

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

Last change on this file since 3398 was 3298, checked in by kanani, 6 years ago
Merge chemistry branch at r3297 to trunk
File size: 42.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |