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, 3 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!
4SUBROUTINE 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!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
302CONTAINS !  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
442TimeLoop: 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
508UntilAccepted: 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
518Stage: 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)*Direction*H
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)/(Direction*H(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) = Direction*H(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!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1252END SUBROUTINE Rosenbrock
1253!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1254
1255
1256
1257
1258
Note: See TracBrowser for help on using the repository browser.