Home


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

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

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

Merge of branch palm4u into trunk
File size: 43.6 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         CALL WCOPY(N,Fcn0,1,Fcn,1)
542      ! istage>1 and a new function evaluation is needed at the current istage
543       ELSEIF ( ros_NewF(istage) ) THEN
544         CALL WCOPY(N,Y,1,Ynew,1)
545         DO j = 1, istage-1
546           CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), &
547            K(N*(j-1)+1:N*j),1,Ynew,1)
548         END DO
549         Tau = T + ros_Alpha(istage)*Direction*H
550         CALL FunTemplate(Tau,Ynew,Fcn)
551         ISTATUS(Nfun) = ISTATUS(Nfun) + 1
552       END IF ! if istage == 1 elseif ros_NewF(istage)
553       CALL WCOPY(N,Fcn,1,K(ioffset+1:ioffset+N),1)
554       DO j = 1, istage-1
555         HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H)
556         CALL WAXPY(N,HC,K(N*(j-1)+1:N*j),1,K(ioffset+1:ioffset+N),1)
557       END DO
558       IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
559         HG = Direction*H*ros_Gamma(istage)
560         CALL WAXPY(N,HG,dFdT,1,K(ioffset+1:ioffset+N),1)
561       END IF
562       CALL ros_Solve(Ghimj, Pivot, K(ioffset+1:ioffset+N))
563
564   END DO Stage
565
566
567!~~~>  Compute the new solution
568   CALL WCOPY(N,Y,1,Ynew,1)
569   DO j=1,ros_S
570         CALL WAXPY(N,ros_M(j),K(N*(j-1)+1:N*j),1,Ynew,1)
571   END DO
572
573!~~~>  Compute the error estimation
574   CALL WSCAL(N,ZERO,Yerr,1)
575   DO j=1,ros_S
576        CALL WAXPY(N,ros_E(j),K(N*(j-1)+1:N*j),1,Yerr,1)
577   END DO
578   Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
579
580!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
581   Fac  = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
582   Hnew = H*Fac
583
584!~~~>  Check the error magnitude and adjust step size
585   ISTATUS(Nstp) = ISTATUS(Nstp) + 1
586   IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN  !~~~> Accept step
587      ISTATUS(Nacc) = ISTATUS(Nacc) + 1
588      CALL WCOPY(N,Ynew,1,Y,1)
589      T = T + Direction*H
590      Hnew = MAX(Hmin,MIN(Hnew,Hmax))
591      IF (RejectLastH) THEN  ! No step size increase after a rejected step
592         Hnew = MIN(Hnew,H)
593      END IF
594      RSTATUS(Nhexit) = H
595      RSTATUS(Nhnew)  = Hnew
596      RSTATUS(Ntexit) = T
597      RejectLastH = .FALSE.
598      RejectMoreH = .FALSE.
599      H = Hnew
600      EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
601   ELSE           !~~~> Reject step
602      IF (RejectMoreH) THEN
603         Hnew = H*FacRej
604      END IF
605      RejectMoreH = RejectLastH
606      RejectLastH = .TRUE.
607      H = Hnew
608      IF (ISTATUS(Nacc) >= 1)  ISTATUS(Nrej) = ISTATUS(Nrej) + 1
609   END IF ! Err <= 1
610
611   END DO UntilAccepted
612
613   END DO TimeLoop
614
615!~~~> Succesful exit
616   IERR = 1  !~~~> The integration was successful
617
618  END SUBROUTINE ros_Integrator
619
620
621!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
622  KPP_REAL FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
623                               AbsTol, RelTol, VectorTol )
624!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
625!~~~> Computes the "scaled norm" of the error vector Yerr
626!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
627   IMPLICIT NONE
628
629! Input arguments
630   KPP_REAL, INTENT(IN) :: Y(N), Ynew(N), &
631          Yerr(N), AbsTol(N), RelTol(N)
632   LOGICAL, INTENT(IN) ::  VectorTol
633! Local variables
634   KPP_REAL :: Err, Scale, Ymax
635   INTEGER  :: i
636   KPP_REAL, PARAMETER :: ZERO = 0.0_dp
637
638   Err = ZERO
639   DO i=1,N
640     Ymax = MAX(ABS(Y(i)),ABS(Ynew(i)))
641     IF (VectorTol) THEN
642       Scale = AbsTol(i)+RelTol(i)*Ymax
643     ELSE
644       Scale = AbsTol(1)+RelTol(1)*Ymax
645     END IF
646     Err = Err+(Yerr(i)/Scale)**2
647   END DO
648   Err  = SQRT(Err/N)
649
650   ros_ErrorNorm = MAX(Err,1.0d-10)
651
652  END FUNCTION ros_ErrorNorm
653
654
655!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
656  SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, &
657                Fcn0, dFdT )
658!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
659!~~~> The time partial derivative of the function by finite differences
660!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
661   IMPLICIT NONE
662
663!~~~> Input arguments
664   KPP_REAL, INTENT(IN) :: T, Roundoff, Y(N), Fcn0(N)
665!~~~> Output arguments
666   KPP_REAL, INTENT(OUT) :: dFdT(N)
667!~~~> Local variables
668   KPP_REAL :: Delta
669   KPP_REAL, PARAMETER :: ONE = 1.0_dp, DeltaMin = 1.0E-6_dp
670
671   Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
672   CALL FunTemplate(T+Delta,Y,dFdT)
673   ISTATUS(Nfun) = ISTATUS(Nfun) + 1
674   CALL WAXPY(N,(-ONE),Fcn0,1,dFdT,1)
675   CALL WSCAL(N,(ONE/Delta),dFdT,1)
676
677  END SUBROUTINE ros_FunTimeDerivative
678
679
680!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
681  SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, &
682             Jac0, Ghimj, Pivot, Singular )
683! --- --- --- --- --- --- --- --- --- --- --- --- ---
684!  Prepares the LHS matrix for stage calculations
685!  1.  Construct Ghimj = 1/(H*ham) - Jac0
686!      "(Gamma H) Inverse Minus Jacobian"
687!  2.  Repeat LU decomposition of Ghimj until successful.
688!       -half the step size if LU decomposition fails and retry
689!       -exit after 5 consecutive fails
690! --- --- --- --- --- --- --- --- --- --- --- --- ---
691   IMPLICIT NONE
692
693!~~~> Input arguments
694#ifdef FULL_ALGEBRA   
695   KPP_REAL, INTENT(IN) ::  Jac0(N,N)
696#else
697   KPP_REAL, INTENT(IN) ::  Jac0(LU_NONZERO)
698#endif   
699   KPP_REAL, INTENT(IN) ::  gam
700   INTEGER, INTENT(IN) ::  Direction
701!~~~> Output arguments
702#ifdef FULL_ALGEBRA   
703   KPP_REAL, INTENT(OUT) :: Ghimj(N,N)
704#else
705   KPP_REAL, INTENT(OUT) :: Ghimj(LU_NONZERO)
706#endif   
707   LOGICAL, INTENT(OUT) ::  Singular
708   INTEGER, INTENT(OUT) ::  Pivot(N)
709!~~~> Inout arguments
710   KPP_REAL, INTENT(INOUT) :: H   ! step size is decreased when LU fails
711!~~~> Local variables
712   INTEGER  :: i, ising, Nconsecutive
713   KPP_REAL :: ghinv
714   KPP_REAL, PARAMETER :: ONE  = 1.0_dp, HALF = 0.5_dp
715
716   Nconsecutive = 0
717   Singular = .TRUE.
718
719   DO WHILE (Singular)
720
721!~~~>    Construct Ghimj = 1/(H*gam) - Jac0
722#ifdef FULL_ALGEBRA   
723     CALL WCOPY(N*N,Jac0,1,Ghimj,1)
724     CALL WSCAL(N*N,(-ONE),Ghimj,1)
725     ghinv = ONE/(Direction*H*gam)
726     DO i=1,N
727       Ghimj(i,i) = Ghimj(i,i)+ghinv
728     END DO
729#else
730     CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1)
731     CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
732     ghinv = ONE/(Direction*H*gam)
733     DO i=1,N
734       Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv
735     END DO
736#endif   
737!~~~>    Compute LU decomposition
738     CALL ros_Decomp( Ghimj, Pivot, ising )
739     IF (ising == 0) THEN
740!~~~>    If successful done
741        Singular = .FALSE.
742     ELSE ! ising .ne. 0
743!~~~>    If unsuccessful half the step size; if 5 consecutive fails then return
744        ISTATUS(Nsng) = ISTATUS(Nsng) + 1
745        Nconsecutive = Nconsecutive+1
746        Singular = .TRUE.
747        PRINT*,'Warning: LU Decomposition returned ising = ',ising
748        IF (Nconsecutive <= 5) THEN ! Less than 5 consecutive failed decompositions
749           H = H*HALF
750        ELSE  ! More than 5 consecutive failed decompositions
751           RETURN
752        END IF  ! Nconsecutive
753      END IF    ! ising
754
755   END DO ! WHILE Singular
756
757  END SUBROUTINE ros_PrepareMatrix
758
759
760!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
761  SUBROUTINE ros_Decomp( A, Pivot, ising )
762!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
763!  Template for the LU decomposition
764!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
765   IMPLICIT NONE
766!~~~> Inout variables
767   KPP_REAL, INTENT(INOUT) :: A(LU_NONZERO)
768!~~~> Output variables
769   INTEGER, INTENT(OUT) :: Pivot(N), ising
770
771#ifdef FULL_ALGEBRA   
772   CALL  DGETRF( N, N, A, N, Pivot, ising )
773#else   
774   CALL KppDecomp ( A, ising )
775   Pivot(1) = 1
776#endif
777   ISTATUS(Ndec) = ISTATUS(Ndec) + 1
778
779  END SUBROUTINE ros_Decomp
780
781
782!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
783  SUBROUTINE ros_Solve( A, Pivot, b )
784!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
785!  Template for the forward/backward substitution (using pre-computed LU decomposition)
786!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
787   IMPLICIT NONE
788!~~~> Input variables
789   KPP_REAL, INTENT(IN) :: A(LU_NONZERO)
790   INTEGER, INTENT(IN) :: Pivot(N)
791!~~~> InOut variables
792   KPP_REAL, INTENT(INOUT) :: b(N)
793
794#ifdef FULL_ALGEBRA   
795   CALL  DGETRS( 'N', N , 1, A, N, Pivot, b, N, 0 )
796#else   
797   CALL KppSolve( A, b )
798#endif
799
800   ISTATUS(Nsol) = ISTATUS(Nsol) + 1
801
802  END SUBROUTINE ros_Solve
803
804
805
806!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
807  SUBROUTINE Ros2
808!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
809! --- AN L-STABLE METHOD, 2 stages, order 2
810!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
811
812   IMPLICIT NONE
813   DOUBLE PRECISION g
814
815    g = 1.0_dp + 1.0_dp/SQRT(2.0_dp)
816    rosMethod = RS2
817!~~~> Name of the method
818    ros_Name = 'ROS-2'
819!~~~> Number of stages
820    ros_S = 2
821
822!~~~> The coefficient matrices A and C are strictly lower triangular.
823!   The lower triangular (subdiagonal) elements are stored in row-wise order:
824!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
825!   The general mapping formula is:
826!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
827!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
828
829    ros_A(1) = (1.0_dp)/g
830    ros_C(1) = (-2.0_dp)/g
831!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
832!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
833    ros_NewF(1) = .TRUE.
834    ros_NewF(2) = .TRUE.
835!~~~> M_i = Coefficients for new step solution
836    ros_M(1)= (3.0_dp)/(2.0_dp*g)
837    ros_M(2)= (1.0_dp)/(2.0_dp*g)
838! E_i = Coefficients for error estimator
839    ros_E(1) = 1.0_dp/(2.0_dp*g)
840    ros_E(2) = 1.0_dp/(2.0_dp*g)
841!~~~> ros_ELO = estimator of local order - the minimum between the
842!    main and the embedded scheme orders plus one
843    ros_ELO = 2.0_dp
844!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
845    ros_Alpha(1) = 0.0_dp
846    ros_Alpha(2) = 1.0_dp
847!~~~> Gamma_i = \sum_j  gamma_{i,j}
848    ros_Gamma(1) = g
849    ros_Gamma(2) =-g
850
851 END SUBROUTINE Ros2
852
853
854!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
855  SUBROUTINE Ros3
856!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
857! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
858!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
859
860   IMPLICIT NONE
861   rosMethod = RS3
862!~~~> Name of the method
863   ros_Name = 'ROS-3'
864!~~~> Number of stages
865   ros_S = 3
866
867!~~~> The coefficient matrices A and C are strictly lower triangular.
868!   The lower triangular (subdiagonal) elements are stored in row-wise order:
869!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
870!   The general mapping formula is:
871!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
872!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
873
874   ros_A(1)= 1.0_dp
875   ros_A(2)= 1.0_dp
876   ros_A(3)= 0.0_dp
877
878   ros_C(1) = -0.10156171083877702091975600115545E+01_dp
879   ros_C(2) =  0.40759956452537699824805835358067E+01_dp
880   ros_C(3) =  0.92076794298330791242156818474003E+01_dp
881!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
882!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
883   ros_NewF(1) = .TRUE.
884   ros_NewF(2) = .TRUE.
885   ros_NewF(3) = .FALSE.
886!~~~> M_i = Coefficients for new step solution
887   ros_M(1) =  0.1E+01_dp
888   ros_M(2) =  0.61697947043828245592553615689730E+01_dp
889   ros_M(3) = -0.42772256543218573326238373806514_dp
890! E_i = Coefficients for error estimator
891   ros_E(1) =  0.5_dp
892   ros_E(2) = -0.29079558716805469821718236208017E+01_dp
893   ros_E(3) =  0.22354069897811569627360909276199_dp
894!~~~> ros_ELO = estimator of local order - the minimum between the
895!    main and the embedded scheme orders plus 1
896   ros_ELO = 3.0_dp
897!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
898   ros_Alpha(1)= 0.0_dp
899   ros_Alpha(2)= 0.43586652150845899941601945119356_dp
900   ros_Alpha(3)= 0.43586652150845899941601945119356_dp
901!~~~> Gamma_i = \sum_j  gamma_{i,j}
902   ros_Gamma(1)= 0.43586652150845899941601945119356_dp
903   ros_Gamma(2)= 0.24291996454816804366592249683314_dp
904   ros_Gamma(3)= 0.21851380027664058511513169485832E+01_dp
905
906  END SUBROUTINE Ros3
907
908!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
909
910
911!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
912  SUBROUTINE Ros4
913!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
914!     L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
915!     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
916!
917!      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
918!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
919!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
920!      SPRINGER-VERLAG (1990)
921!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
922
923   IMPLICIT NONE
924
925   rosMethod = RS4
926!~~~> Name of the method
927   ros_Name = 'ROS-4'
928!~~~> Number of stages
929   ros_S = 4
930
931!~~~> The coefficient matrices A and C are strictly lower triangular.
932!   The lower triangular (subdiagonal) elements are stored in row-wise order:
933!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
934!   The general mapping formula is:
935!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
936!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
937
938   ros_A(1) = 0.2000000000000000E+01_dp
939   ros_A(2) = 0.1867943637803922E+01_dp
940   ros_A(3) = 0.2344449711399156_dp
941   ros_A(4) = ros_A(2)
942   ros_A(5) = ros_A(3)
943   ros_A(6) = 0.0_dp
944
945   ros_C(1) =-0.7137615036412310E+01_dp
946   ros_C(2) = 0.2580708087951457E+01_dp
947   ros_C(3) = 0.6515950076447975_dp
948   ros_C(4) =-0.2137148994382534E+01_dp
949   ros_C(5) =-0.3214669691237626_dp
950   ros_C(6) =-0.6949742501781779_dp
951!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
952!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
953   ros_NewF(1)  = .TRUE.
954   ros_NewF(2)  = .TRUE.
955   ros_NewF(3)  = .TRUE.
956   ros_NewF(4)  = .FALSE.
957!~~~> M_i = Coefficients for new step solution
958   ros_M(1) = 0.2255570073418735E+01_dp
959   ros_M(2) = 0.2870493262186792_dp
960   ros_M(3) = 0.4353179431840180_dp
961   ros_M(4) = 0.1093502252409163E+01_dp
962!~~~> E_i  = Coefficients for error estimator
963   ros_E(1) =-0.2815431932141155_dp
964   ros_E(2) =-0.7276199124938920E-01_dp
965   ros_E(3) =-0.1082196201495311_dp
966   ros_E(4) =-0.1093502252409163E+01_dp
967!~~~> ros_ELO  = estimator of local order - the minimum between the
968!    main and the embedded scheme orders plus 1
969   ros_ELO  = 4.0_dp
970!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
971   ros_Alpha(1) = 0.0_dp
972   ros_Alpha(2) = 0.1145640000000000E+01_dp
973   ros_Alpha(3) = 0.6552168638155900_dp
974   ros_Alpha(4) = ros_Alpha(3)
975!~~~> Gamma_i = \sum_j  gamma_{i,j}
976   ros_Gamma(1) = 0.5728200000000000_dp
977   ros_Gamma(2) =-0.1769193891319233E+01_dp
978   ros_Gamma(3) = 0.7592633437920482_dp
979   ros_Gamma(4) =-0.1049021087100450_dp
980
981  END SUBROUTINE Ros4
982
983!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
984  SUBROUTINE Rodas3
985!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
986! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
987!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
988
989   IMPLICIT NONE
990
991   rosMethod = RD3
992!~~~> Name of the method
993   ros_Name = 'RODAS-3'
994!~~~> Number of stages
995   ros_S = 4
996
997!~~~> The coefficient matrices A and C are strictly lower triangular.
998!   The lower triangular (subdiagonal) elements are stored in row-wise order:
999!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1000!   The general mapping formula is:
1001!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1002!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1003
1004   ros_A(1) = 0.0_dp
1005   ros_A(2) = 2.0_dp
1006   ros_A(3) = 0.0_dp
1007   ros_A(4) = 2.0_dp
1008   ros_A(5) = 0.0_dp
1009   ros_A(6) = 1.0_dp
1010
1011   ros_C(1) = 4.0_dp
1012   ros_C(2) = 1.0_dp
1013   ros_C(3) =-1.0_dp
1014   ros_C(4) = 1.0_dp
1015   ros_C(5) =-1.0_dp
1016   ros_C(6) =-(8.0_dp/3.0_dp)
1017
1018!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1019!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1020   ros_NewF(1)  = .TRUE.
1021   ros_NewF(2)  = .FALSE.
1022   ros_NewF(3)  = .TRUE.
1023   ros_NewF(4)  = .TRUE.
1024!~~~> M_i = Coefficients for new step solution
1025   ros_M(1) = 2.0_dp
1026   ros_M(2) = 0.0_dp
1027   ros_M(3) = 1.0_dp
1028   ros_M(4) = 1.0_dp
1029!~~~> E_i  = Coefficients for error estimator
1030   ros_E(1) = 0.0_dp
1031   ros_E(2) = 0.0_dp
1032   ros_E(3) = 0.0_dp
1033   ros_E(4) = 1.0_dp
1034!~~~> ros_ELO  = estimator of local order - the minimum between the
1035!    main and the embedded scheme orders plus 1
1036   ros_ELO  = 3.0_dp
1037!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1038   ros_Alpha(1) = 0.0_dp
1039   ros_Alpha(2) = 0.0_dp
1040   ros_Alpha(3) = 1.0_dp
1041   ros_Alpha(4) = 1.0_dp
1042!~~~> Gamma_i = \sum_j  gamma_{i,j}
1043   ros_Gamma(1) = 0.5_dp
1044   ros_Gamma(2) = 1.5_dp
1045   ros_Gamma(3) = 0.0_dp
1046   ros_Gamma(4) = 0.0_dp
1047
1048  END SUBROUTINE Rodas3
1049
1050!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1051  SUBROUTINE Rodas4
1052!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1053!     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
1054!
1055!      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
1056!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1057!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1058!      SPRINGER-VERLAG (1996)
1059!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1060
1061   IMPLICIT NONE
1062
1063    rosMethod = RD4
1064!~~~> Name of the method
1065    ros_Name = 'RODAS-4'
1066!~~~> Number of stages
1067    ros_S = 6
1068
1069!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1070    ros_Alpha(1) = 0.000_dp
1071    ros_Alpha(2) = 0.386_dp
1072    ros_Alpha(3) = 0.210_dp
1073    ros_Alpha(4) = 0.630_dp
1074    ros_Alpha(5) = 1.000_dp
1075    ros_Alpha(6) = 1.000_dp
1076
1077!~~~> Gamma_i = \sum_j  gamma_{i,j}
1078    ros_Gamma(1) = 0.2500000000000000_dp
1079    ros_Gamma(2) =-0.1043000000000000_dp
1080    ros_Gamma(3) = 0.1035000000000000_dp
1081    ros_Gamma(4) =-0.3620000000000023E-01_dp
1082    ros_Gamma(5) = 0.0_dp
1083    ros_Gamma(6) = 0.0_dp
1084
1085!~~~> The coefficient matrices A and C are strictly lower triangular.
1086!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1087!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1088!   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1089!                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1090
1091    ros_A(1) = 0.1544000000000000E+01_dp
1092    ros_A(2) = 0.9466785280815826_dp
1093    ros_A(3) = 0.2557011698983284_dp
1094    ros_A(4) = 0.3314825187068521E+01_dp
1095    ros_A(5) = 0.2896124015972201E+01_dp
1096    ros_A(6) = 0.9986419139977817_dp
1097    ros_A(7) = 0.1221224509226641E+01_dp
1098    ros_A(8) = 0.6019134481288629E+01_dp
1099    ros_A(9) = 0.1253708332932087E+02_dp
1100    ros_A(10) =-0.6878860361058950_dp
1101    ros_A(11) = ros_A(7)
1102    ros_A(12) = ros_A(8)
1103    ros_A(13) = ros_A(9)
1104    ros_A(14) = ros_A(10)
1105    ros_A(15) = 1.0_dp
1106
1107    ros_C(1) =-0.5668800000000000E+01_dp
1108    ros_C(2) =-0.2430093356833875E+01_dp
1109    ros_C(3) =-0.2063599157091915_dp
1110    ros_C(4) =-0.1073529058151375_dp
1111    ros_C(5) =-0.9594562251023355E+01_dp
1112    ros_C(6) =-0.2047028614809616E+02_dp
1113    ros_C(7) = 0.7496443313967647E+01_dp
1114    ros_C(8) =-0.1024680431464352E+02_dp
1115    ros_C(9) =-0.3399990352819905E+02_dp
1116    ros_C(10) = 0.1170890893206160E+02_dp
1117    ros_C(11) = 0.8083246795921522E+01_dp
1118    ros_C(12) =-0.7981132988064893E+01_dp
1119    ros_C(13) =-0.3152159432874371E+02_dp
1120    ros_C(14) = 0.1631930543123136E+02_dp
1121    ros_C(15) =-0.6058818238834054E+01_dp
1122
1123!~~~> M_i = Coefficients for new step solution
1124    ros_M(1) = ros_A(7)
1125    ros_M(2) = ros_A(8)
1126    ros_M(3) = ros_A(9)
1127    ros_M(4) = ros_A(10)
1128    ros_M(5) = 1.0_dp
1129    ros_M(6) = 1.0_dp
1130
1131!~~~> E_i  = Coefficients for error estimator
1132    ros_E(1) = 0.0_dp
1133    ros_E(2) = 0.0_dp
1134    ros_E(3) = 0.0_dp
1135    ros_E(4) = 0.0_dp
1136    ros_E(5) = 0.0_dp
1137    ros_E(6) = 1.0_dp
1138
1139!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1140!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1141    ros_NewF(1) = .TRUE.
1142    ros_NewF(2) = .TRUE.
1143    ros_NewF(3) = .TRUE.
1144    ros_NewF(4) = .TRUE.
1145    ros_NewF(5) = .TRUE.
1146    ros_NewF(6) = .TRUE.
1147
1148!~~~> ros_ELO  = estimator of local order - the minimum between the
1149!        main and the embedded scheme orders plus 1
1150    ros_ELO = 4.0_dp
1151
1152  END SUBROUTINE Rodas4
1153!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1154  SUBROUTINE Rang3
1155!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1156! STIFFLY-STABLE W METHOD OF ORDER 3, WITH 4 STAGES
1157!
1158! J. RANG and L. ANGERMANN
1159! NEW ROSENBROCK W-METHODS OF ORDER 3
1160! FOR PARTIAL DIFFERENTIAL ALGEBRAIC
1161!        EQUATIONS OF INDEX 1                                             
1162! BIT Numerical Mathematics (2005) 45: 761-787
1163!  DOI: 10.1007/s10543-005-0035-y
1164! Table 4.1-4.2
1165!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1166
1167   IMPLICIT NONE
1168
1169    rosMethod = RG3
1170!~~~> Name of the method
1171    ros_Name = 'RANG-3'
1172!~~~> Number of stages
1173    ros_S = 4
1174
1175    ros_A(1) = 5.09052051067020d+00;
1176    ros_A(2) = 5.09052051067020d+00;
1177    ros_A(3) = 0.0d0;
1178    ros_A(4) = 4.97628111010787d+00;
1179    ros_A(5) = 2.77268164715849d-02;
1180    ros_A(6) = 2.29428036027904d-01;
1181
1182    ros_C(1) = -1.16790812312283d+01;
1183    ros_C(2) = -1.64057326467367d+01;
1184    ros_C(3) = -2.77268164715850d-01;
1185    ros_C(4) = -8.38103960500476d+00;
1186    ros_C(5) = -8.48328409199343d-01;
1187    ros_C(6) =  2.87009860433106d-01;
1188
1189    ros_M(1) =  5.22582761233094d+00;
1190    ros_M(2) = -5.56971148154165d-01;
1191    ros_M(3) =  3.57979469353645d-01;
1192    ros_M(4) =  1.72337398521064d+00;
1193
1194    ros_E(1) = -5.16845212784040d+00;
1195    ros_E(2) = -1.26351942603842d+00;
1196    ros_E(3) = -1.11022302462516d-16;
1197    ros_E(4) =  2.22044604925031d-16;
1198
1199    ros_Alpha(1) = 0.0d00;
1200    ros_Alpha(2) = 2.21878746765329d+00;
1201    ros_Alpha(3) = 2.21878746765329d+00;
1202    ros_Alpha(4) = 1.55392337535788d+00;
1203
1204    ros_Gamma(1) =  4.35866521508459d-01;
1205    ros_Gamma(2) = -1.78292094614483d+00;
1206    ros_Gamma(3) = -2.46541900496934d+00;
1207    ros_Gamma(4) = -8.05529997906370d-01;
1208
1209
1210!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1211!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1212    ros_NewF(1) = .TRUE.
1213    ros_NewF(2) = .TRUE.
1214    ros_NewF(3) = .TRUE.
1215    ros_NewF(4) = .TRUE.
1216
1217!~~~> ros_ELO  = estimator of local order - the minimum between the
1218!        main and the embedded scheme orders plus 1
1219    ros_ELO = 3.0_dp
1220
1221  END SUBROUTINE Rang3
1222!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1223
1224!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1225!   End of the set of internal Rosenbrock subroutines
1226!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1227END SUBROUTINE Rosenbrock
1228!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1229
1230
1231!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1232SUBROUTINE FunTemplate( T, Y, Ydot )
1233!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1234!  Template for the ODE function call.
1235!  Updates the rate coefficients (and possibly the fixed species) at each call
1236!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1237 USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO
1238 USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME
1239 USE KPP_ROOT_Function, ONLY: Fun
1240 USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST
1241!~~~> Input variables
1242   KPP_REAL :: T, Y(NVAR)
1243!~~~> Output variables
1244   KPP_REAL :: Ydot(NVAR)
1245!~~~> Local variables
1246   KPP_REAL :: Told
1247
1248   Told = TIME
1249   TIME = T
1250   CALL Update_SUN()
1251   CALL Update_RCONST()
1252   CALL Fun( Y, FIX, RCONST, Ydot )
1253   TIME = Told
1254
1255END SUBROUTINE FunTemplate
1256
1257
1258!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1259SUBROUTINE JacTemplate( T, Y, Jcb )
1260!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1261!  Template for the ODE Jacobian call.
1262!  Updates the rate coefficients (and possibly the fixed species) at each call
1263!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1264 USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO
1265 USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME
1266 USE KPP_ROOT_Jacobian, ONLY: Jac_SP, LU_IROW, LU_ICOL
1267 USE KPP_ROOT_LinearAlgebra
1268 USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST
1269!~~~> Input variables
1270    KPP_REAL :: T, Y(NVAR)
1271!~~~> Output variables
1272#ifdef FULL_ALGEBRA   
1273    KPP_REAL :: JV(LU_NONZERO), Jcb(NVAR,NVAR)
1274#else
1275    KPP_REAL :: Jcb(LU_NONZERO)
1276#endif   
1277!~~~> Local variables
1278    KPP_REAL :: Told
1279#ifdef FULL_ALGEBRA   
1280    INTEGER :: i, j
1281#endif   
1282
1283    Told = TIME
1284    TIME = T
1285    CALL Update_SUN()
1286    CALL Update_RCONST()
1287#ifdef FULL_ALGEBRA   
1288    CALL Jac_SP(Y, FIX, RCONST, JV)
1289    DO j=1,NVAR
1290      DO i=1,NVAR
1291         Jcb(i,j) = 0.0_dp
1292      END DO
1293    END DO
1294    DO i=1,LU_NONZERO
1295       Jcb(LU_IROW(i),LU_ICOL(i)) = JV(i)
1296    END DO
1297#else
1298    CALL Jac_SP( Y, FIX, RCONST, Jcb )
1299#endif   
1300    TIME = Told
1301
1302END SUBROUTINE JacTemplate
1303
1304END MODULE KPP_ROOT_Integrator
1305
1306
1307
1308
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |