Home

Context Navigation

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

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

Rev	Line
[3185]	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)
[3260]	167
[3185]	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 |