Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/rosenbrock.f90 @ 3780

Last change on this file since 3780 was 2696, checked in by kanani, 7 years ago

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