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

Last change on this file since 3458 was 3458, checked in by kanani, 6 years ago

Reintegrated fixes/changes from branch chemistry

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