source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/sdirk.f90 @ 4794

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

Merge of branch palm4u into trunk

File size: 42.4 KB
Line 
1!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
2!  SDIRK - Singly-Diagonally-Implicit Runge-Kutta methods                 !
3!            * Sdirk 2a, 2b: L-stable, 2 stages, order 2                  !
4!            * Sdirk 3a:     L-stable, 3 stages, order 2, adj-invariant   !
5!            * Sdirk 4a, 4b: L-stable, 5 stages, order 4                  !
6!  By default the code employs the KPP sparse linear algebra routines     !
7!  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
8!                                                                         !
9!    (C)  Adrian Sandu, July 2005                                         !
10!    Virginia Polytechnic Institute and State University                  !
11!    Contact: sandu@cs.vt.edu                                             !
12!    Revised by Philipp Miehe and Adrian Sandu, May 2006                  !
13!    This implementation is part of KPP - the Kinetic PreProcessor        !
14!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
15
16MODULE KPP_ROOT_Integrator
17
18  USE KPP_ROOT_Precision
19  USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME
20  USE KPP_ROOT_Parameters, ONLY: NVAR, NSPEC, NFIX, LU_NONZERO
21  USE KPP_ROOT_JacobianSP, ONLY: LU_DIAG
22  USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp,    &
23               KppSolve, Set2zero, WLAMCH, WCOPY, WAXPY, WSCAL, WADD
24 
25  IMPLICIT NONE
26  PUBLIC
27  SAVE
28 
29!~~~>  Statistics on the work performed by the SDIRK method
30  INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4,  &
31           Nrej=5, Ndec=6, Nsol=7, Nsng=8,               &
32           Ntexit=1, Nhexit=2, Nhnew=3
33                 
34CONTAINS
35
36SUBROUTINE INTEGRATE( TIN, TOUT, &
37  ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, Ierr_U )
38
39   USE KPP_ROOT_Parameters
40   USE KPP_ROOT_Global
41   IMPLICIT NONE
42
43   KPP_REAL, INTENT(IN) :: TIN  ! Start Time
44   KPP_REAL, INTENT(IN) :: TOUT ! End Time
45   ! Optional input parameters and statistics
46   INTEGER,       INTENT(IN),  OPTIONAL :: ICNTRL_U(20)
47   KPP_REAL, INTENT(IN),  OPTIONAL :: RCNTRL_U(20)
48   INTEGER,       INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
49   KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
50   INTEGER,       INTENT(OUT), OPTIONAL :: Ierr_U
51
52   INTEGER, SAVE :: Ntotal = 0
53   KPP_REAL :: RCNTRL(20), RSTATUS(20), T1, T2
54   INTEGER       :: ICNTRL(20), ISTATUS(20), Ierr
55
56   ICNTRL(:)  = 0
57   RCNTRL(:)  = 0.0_dp
58   ISTATUS(:) = 0
59   RSTATUS(:) = 0.0_dp
60
61   !~~~> fine-tune the integrator:
62   ICNTRL(2) = 0        ! 0 - vector tolerances, 1 - scalar tolerances
63   ICNTRL(6) = 0        ! starting values of Newton iterations: interpolated (0), zero (1)
64  ! If optional parameters are given, and if they are >0,
65   ! then they overwrite default settings.
66   IF (PRESENT(ICNTRL_U)) THEN
67     WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)
68   END IF
69   IF (PRESENT(RCNTRL_U)) THEN
70     WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:)
71   END IF
72   
73   
74   T1 = TIN; T2 = TOUT
75   CALL SDIRK( NVAR,T1,T2,VAR,RTOL,ATOL,          &
76               RCNTRL,ICNTRL,RSTATUS,ISTATUS,Ierr )
77
78   !~~~> Debug option: print number of steps
79   ! Ntotal = Ntotal + ISTATUS(Nstp)
80   ! PRINT*,'NSTEPS=',ISTATUS(Nstp), '(',Ntotal,')','   O3=',VAR(ind_O3), &
81   !       '  NO2=',VAR(ind_NO2)
82
83   IF (Ierr < 0) THEN
84        PRINT *,'SDIRK: Unsuccessful exit at T=',TIN,' (Ierr=',Ierr,')'
85   ENDIF
86   
87   ! if optional parameters are given for output they to return information
88   IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:)
89   IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(:)
90   IF (PRESENT(Ierr_U))    Ierr_U       = Ierr
91
92   END SUBROUTINE INTEGRATE
93
94
95!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
96   SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol,     &
97                       RCNTRL, ICNTRL, RSTATUS, ISTATUS, Ierr)
98!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
99!
100!    Solves the system y'=F(t,y) using a Singly-Diagonally-Implicit
101!    Runge-Kutta (SDIRK) method.
102!
103!    This implementation is based on the book and the code Sdirk4:
104!
105!      E. Hairer and G. Wanner
106!      "Solving ODEs II. Stiff and differential-algebraic problems".
107!      Springer series in computational mathematics, Springer-Verlag, 1996.
108!    This code is based on the SDIRK4 routine in the above book.
109!
110!    Methods:
111!            * Sdirk 2a, 2b: L-stable, 2 stages, order 2                 
112!            * Sdirk 3a:     L-stable, 3 stages, order 2, adjoint-invariant   
113!            * Sdirk 4a, 4b: L-stable, 5 stages, order 4                 
114!
115!    (C)  Adrian Sandu, July 2005
116!    Virginia Polytechnic Institute and State University
117!    Contact: sandu@cs.vt.edu
118!    Revised by Philipp Miehe and Adrian Sandu, May 2006                 
119!    This implementation is part of KPP - the Kinetic PreProcessor
120!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
121!
122!~~~>   INPUT ARGUMENTS:
123!
124!-     Y(NVAR)    = vector of initial conditions (at T=Tinitial)
125!-    [Tinitial,Tfinal]  = time range of integration
126!     (if Tinitial>Tfinal the integration is performed backwards in time)
127!-    RelTol, AbsTol = user precribed accuracy
128!- SUBROUTINE ode_Fun( T, Y, Ydot ) = ODE function,
129!                       returns Ydot = Y' = F(T,Y)
130!- SUBROUTINE ode_Fun( T, Y, Ydot ) = Jacobian of the ODE function,
131!                       returns Jcb = dF/dY
132!-    ICNTRL(1:20)    = integer inputs parameters
133!-    RCNTRL(1:20)    = real inputs parameters
134!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
135!
136!~~~>     OUTPUT ARGUMENTS:
137!
138!-    Y(NVAR)         -> vector of final states (at T->Tfinal)
139!-    ISTATUS(1:20)   -> integer output parameters
140!-    RSTATUS(1:20)   -> real output parameters
141!-    Ierr            -> job status upon return
142!                        success (positive value) or
143!                        failure (negative value)
144!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
145!
146!~~~>     INPUT PARAMETERS:
147!
148!    Note: For input parameters equal to zero the default values of the
149!       corresponding variables are used.
150!
151!    Note: For input parameters equal to zero the default values of the
152!          corresponding variables are used.
153!~~~> 
154!    ICNTRL(1) = not used
155!
156!    ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
157!              = 1: AbsTol, RelTol are scalars
158!
159!    ICNTRL(3) = Method
160!
161!    ICNTRL(4)  -> maximum number of integration steps
162!        For ICNTRL(4)=0 the default value of 100000 is used
163!
164!    ICNTRL(5)  -> maximum number of Newton iterations
165!        For ICNTRL(5)=0 the default value of 8 is used
166!
167!    ICNTRL(6)  -> starting values of Newton iterations:
168!        ICNTRL(6)=0 : starting values are interpolated (the default)
169!        ICNTRL(6)=1 : starting values are zero
170!
171!~~~>  Real parameters
172!
173!    RCNTRL(1)  -> Hmin, lower bound for the integration step size
174!                  It is strongly recommended to keep Hmin = ZERO
175!    RCNTRL(2)  -> Hmax, upper bound for the integration step size
176!    RCNTRL(3)  -> Hstart, starting value for the integration step size
177!
178!    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
179!    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
180!    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
181!                 (default=0.1)
182!    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
183!                  than the predicted value  (default=0.9)
184!    RCNTRL(8)  -> ThetaMin. If Newton convergence rate smaller
185!                  than ThetaMin the Jacobian is not recomputed;
186!                  (default=0.001)
187!    RCNTRL(9)  -> NewtonTol, stopping criterion for Newton's method
188!                  (default=0.03)
189!    RCNTRL(10) -> Qmin
190!    RCNTRL(11) -> Qmax. If Qmin < Hnew/Hold < Qmax, then the
191!                  step size is kept constant and the LU factorization
192!                  reused (default Qmin=1, Qmax=1.2)
193!
194!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
195!
196!~~~>     OUTPUT PARAMETERS:
197!
198!    Note: each call to Rosenbrock adds the current no. of fcn calls
199!      to previous value of ISTATUS(1), and similar for the other params.
200!      Set ISTATUS(1:10) = 0 before call to avoid this accumulation.
201!
202!    ISTATUS(1) = No. of function calls
203!    ISTATUS(2) = No. of jacobian calls
204!    ISTATUS(3) = No. of steps
205!    ISTATUS(4) = No. of accepted steps
206!    ISTATUS(5) = No. of rejected steps (except at the beginning)
207!    ISTATUS(6) = No. of LU decompositions
208!    ISTATUS(7) = No. of forward/backward substitutions
209!    ISTATUS(8) = No. of singular matrix decompositions
210!
211!    RSTATUS(1)  -> Texit, the time corresponding to the
212!                   computed Y upon return
213!    RSTATUS(2)  -> Hexit,last accepted step before return
214!    RSTATUS(3)  -> Hnew, last predicted step before return
215!        For multiple restarts, use Hnew as Hstart in the following run
216!
217!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
218      IMPLICIT NONE
219
220! Arguments     
221      INTEGER, INTENT(IN)          :: N, ICNTRL(20)
222      KPP_REAL, INTENT(IN)    :: Tinitial, Tfinal, &
223                    RelTol(NVAR), AbsTol(NVAR), RCNTRL(20)
224      KPP_REAL, INTENT(INOUT) :: Y(NVAR)
225      INTEGER, INTENT(OUT)         :: Ierr
226      INTEGER, INTENT(INOUT)       :: ISTATUS(20) 
227      KPP_REAL, INTENT(OUT)   :: RSTATUS(20)
228       
229!~~~>  SDIRK method coefficients, up to 5 stages
230      INTEGER, PARAMETER :: Smax = 5
231      INTEGER, PARAMETER :: S2A=1, S2B=2, S3A=3, S4A=4, S4B=5
232      KPP_REAL :: rkGamma, rkA(Smax,Smax), rkB(Smax), rkC(Smax), &
233                       rkD(Smax),  rkE(Smax), rkBhat(Smax), rkELO,    &
234                       rkAlpha(Smax,Smax), rkTheta(Smax,Smax)
235      INTEGER :: sdMethod, rkS ! The number of stages
236! Local variables     
237      LOGICAL :: StartNewton
238      KPP_REAL :: Hmin, Hmax, Hstart, Roundoff,    &
239                       FacMin, Facmax, FacSafe, FacRej, &
240                       ThetaMin, NewtonTol, Qmin, Qmax
241      INTEGER       :: ITOL, NewtonMaxit, Max_no_steps, i
242      KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
243
244       
245      ISTATUS(1:20) = 0
246      RSTATUS(1:20) = ZERO
247      Ierr = 0
248
249!~~~>  For Scalar tolerances (ICNTRL(2).NE.0)  the code uses AbsTol(1) and RelTol(1)
250!   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
251      IF (ICNTRL(2) == 0) THEN
252         ITOL = 1
253      ELSE
254         ITOL = 0
255      END IF
256
257!~~~> ICNTRL(3) - method selection       
258      SELECT CASE (ICNTRL(3))
259      CASE (0,1)
260         CALL Sdirk2a
261      CASE (2)
262         CALL Sdirk2b
263      CASE (3)
264         CALL Sdirk3a
265      CASE (4)
266         CALL Sdirk4a
267      CASE (5)
268         CALL Sdirk4b
269      CASE DEFAULT
270         CALL Sdirk2a
271      END SELECT
272     
273!~~~>   The maximum number of time steps admitted
274      IF (ICNTRL(4) == 0) THEN
275         Max_no_steps = 200000
276      ELSEIF (ICNTRL(4) > 0) THEN
277         Max_no_steps=ICNTRL(4)
278      ELSE
279         PRINT * ,'User-selected ICNTRL(4)=',ICNTRL(4)
280         CALL SDIRK_ErrorMsg(-1,Tinitial,ZERO,Ierr)
281   END IF
282
283!~~~> The maximum number of Newton iterations admitted
284      IF(ICNTRL(5) == 0)THEN
285         NewtonMaxit=8
286      ELSE
287         NewtonMaxit=ICNTRL(5)
288         IF(NewtonMaxit <= 0)THEN
289             PRINT * ,'User-selected ICNTRL(5)=',ICNTRL(5)
290             CALL SDIRK_ErrorMsg(-2,Tinitial,ZERO,Ierr)
291         END IF
292      END IF
293
294!~~~> StartNewton:  Use extrapolation for starting values of Newton iterations
295      IF (ICNTRL(6) == 0) THEN
296         StartNewton = .TRUE.
297      ELSE
298         StartNewton = .FALSE.
299      END IF     
300
301!~~~>  Unit roundoff (1+Roundoff>1)
302      Roundoff = WLAMCH('E')
303
304!~~~>  Lower bound on the step size: (positive value)
305      IF (RCNTRL(1) == ZERO) THEN
306         Hmin = ZERO
307      ELSEIF (RCNTRL(1) > ZERO) THEN
308         Hmin = RCNTRL(1)
309      ELSE
310         PRINT * , 'User-selected RCNTRL(1)=', RCNTRL(1)
311         CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr)
312      END IF
313   
314!~~~>  Upper bound on the step size: (positive value)
315      IF (RCNTRL(2) == ZERO) THEN
316         Hmax = ABS(Tfinal-Tinitial)
317      ELSEIF (RCNTRL(2) > ZERO) THEN
318         Hmax = MIN(ABS(RCNTRL(2)),ABS(Tfinal-Tinitial))
319      ELSE
320         PRINT * , 'User-selected RCNTRL(2)=', RCNTRL(2)
321         CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr)
322      END IF
323   
324!~~~>  Starting step size: (positive value)
325      IF (RCNTRL(3) == ZERO) THEN
326         Hstart = MAX(Hmin,Roundoff)
327      ELSEIF (RCNTRL(3) > ZERO) THEN
328         Hstart = MIN(ABS(RCNTRL(3)),ABS(Tfinal-Tinitial))
329      ELSE
330         PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
331         CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr)
332      END IF
333   
334!~~~>  Step size can be changed s.t.  FacMin < Hnew/Hexit < FacMax
335      IF (RCNTRL(4) == ZERO) THEN
336         FacMin = 0.2_dp
337      ELSEIF (RCNTRL(4) > ZERO) THEN
338         FacMin = RCNTRL(4)
339      ELSE
340         PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
341         CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr)
342      END IF
343      IF (RCNTRL(5) == ZERO) THEN
344         FacMax = 10.0_dp
345      ELSEIF (RCNTRL(5) > ZERO) THEN
346         FacMax = RCNTRL(5)
347      ELSE
348         PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
349         CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr)
350      END IF
351!~~~>   FacRej: Factor to decrease step after 2 succesive rejections
352      IF (RCNTRL(6) == ZERO) THEN
353         FacRej = 0.1_dp
354      ELSEIF (RCNTRL(6) > ZERO) THEN
355         FacRej = RCNTRL(6)
356      ELSE
357         PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
358         CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr)
359      END IF
360!~~~>   FacSafe: Safety Factor in the computation of new step size
361      IF (RCNTRL(7) == ZERO) THEN
362         FacSafe = 0.9_dp
363      ELSEIF (RCNTRL(7) > ZERO) THEN
364         FacSafe = RCNTRL(7)
365      ELSE
366         PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
367         CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr)
368      END IF
369
370!~~~> ThetaMin: decides whether the Jacobian should be recomputed
371      IF(RCNTRL(8) == 0.D0)THEN
372         ThetaMin = 1.0d-3
373      ELSE
374         ThetaMin = RCNTRL(8)
375      END IF
376
377!~~~> Stopping criterion for Newton's method
378      IF(RCNTRL(9) == ZERO)THEN
379         NewtonTol = 3.0d-2
380      ELSE
381         NewtonTol = RCNTRL(9)
382      END IF
383
384!~~~> Qmin, Qmax: IF Qmin < Hnew/Hold < Qmax, STEP SIZE = CONST.
385      IF(RCNTRL(10) == ZERO)THEN
386         Qmin=ONE
387      ELSE
388         Qmin=RCNTRL(10)
389      END IF
390      IF(RCNTRL(11) == ZERO)THEN
391         Qmax=1.2D0
392      ELSE
393         Qmax=RCNTRL(11)
394      END IF
395
396!~~~>  Check if tolerances are reasonable
397      IF (ITOL == 0) THEN
398         IF (AbsTol(1) <= ZERO .OR. RelTol(1) <= 10.D0*Roundoff) THEN
399            PRINT * , ' Scalar AbsTol = ',AbsTol(1)
400            PRINT * , ' Scalar RelTol = ',RelTol(1)
401            CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,Ierr)
402         END IF
403      ELSE
404         DO i=1,N
405            IF (AbsTol(i) <= 0.D0.OR.RelTol(i) <= 10.D0*Roundoff) THEN
406              PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
407              PRINT * , ' RelTol(',i,') = ',RelTol(i)
408              CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,Ierr)
409            END IF
410         END DO
411      END IF
412   
413    IF (Ierr < 0) RETURN
414
415    CALL SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr )
416
417
418   !END SUBROUTINE SDIRK
419
420!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
421   CONTAINS  !  PROCEDURES INTERNAL TO SDIRK
422!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
423
424
425!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
426   SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr )
427!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
428
429      USE KPP_ROOT_Parameters
430      IMPLICIT NONE
431
432!~~~> Arguments:     
433      INTEGER, INTENT(IN) :: N
434      KPP_REAL, INTENT(INOUT) :: Y(NVAR)
435      KPP_REAL, INTENT(IN) :: Tinitial, Tfinal
436      INTEGER, INTENT(OUT) :: Ierr
437     
438!~~~> Local variables:     
439      KPP_REAL :: Z(NVAR,Smax), G(NVAR), TMP(NVAR),        &   
440                       NewtonRate, SCAL(NVAR), RHS(NVAR),       &
441                       T, H, Theta, Hratio, NewtonPredictedErr, &
442                       Qnewton, Err, Fac, Hnew,Tdirection,      &
443                       NewtonIncrement, NewtonIncrementOld
444      INTEGER :: j, IER, istage, NewtonIter, IP(NVAR)
445      LOGICAL :: Reject, FirstStep, SkipJac, SkipLU, NewtonDone
446     
447#ifdef FULL_ALGEBRA     
448      KPP_REAL FJAC(NVAR,NVAR),  E(NVAR,NVAR)
449#else     
450      KPP_REAL FJAC(LU_NONZERO), E(LU_NONZERO)
451#endif     
452      KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
453
454
455!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
456!~~~>  Initializations
457!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
458
459      T = Tinitial
460      Tdirection = SIGN(ONE,Tfinal-Tinitial)
461      H = MAX(ABS(Hmin),ABS(Hstart))
462      IF (ABS(H) <= 10.D0*Roundoff) H=1.0D-6
463      H=MIN(ABS(H),Hmax)
464      H=SIGN(H,Tdirection)
465      SkipLU =.FALSE.
466      SkipJac = .FALSE.
467      Reject=.FALSE.
468      FirstStep=.TRUE.
469
470      CALL SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL)
471
472!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
473!~~~>  Time loop begins
474!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
475Tloop: DO WHILE ( (Tfinal-T)*Tdirection - Roundoff > ZERO )
476
477
478!~~~>  Compute E = 1/(h*gamma)-Jac and its LU decomposition
479      IF ( .NOT.SkipLU ) THEN ! This time around skip the Jac update and LU
480         CALL SDIRK_PrepareMatrix ( H, T, Y, FJAC, &
481                   SkipJac, SkipLU, E, IP, Reject, IER )
482         IF (IER /= 0) THEN
483             CALL SDIRK_ErrorMsg(-8,T,H,Ierr); RETURN
484         END IF
485      END IF     
486
487      IF (ISTATUS(Nstp) > Max_no_steps) THEN
488             CALL SDIRK_ErrorMsg(-6,T,H,Ierr); RETURN
489      END IF   
490      IF ( (T+0.1d0*H == T) .OR. (ABS(H) <= Roundoff) ) THEN
491             CALL SDIRK_ErrorMsg(-7,T,H,Ierr); RETURN
492      END IF   
493
494stages:DO istage = 1, rkS
495
496!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
497!~~~>  Simplified Newton iterations
498!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
499
500!~~~>  Starting values for Newton iterations
501       CALL Set2zero(N,Z(1,istage))
502       
503!~~~>   Prepare the loop-independent part of the right-hand side
504       CALL Set2zero(N,G)
505       IF (istage > 1) THEN
506           DO j = 1, istage-1
507               ! Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj)
508               CALL WAXPY(N,rkTheta(istage,j),Z(1,j),1,G,1)
509               ! Zi(:) = sum_j Alpha(i,j)*Zj(:)
510               IF (StartNewton) THEN
511                  CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1)
512               END IF
513           END DO
514       END IF
515
516       !~~~>  Initializations for Newton iteration
517       NewtonDone = .FALSE.
518       Fac = 0.5d0 ! Step reduction factor if too many iterations
519           
520NewtonLoop:DO NewtonIter = 1, NewtonMaxit
521
522!~~~>   Prepare the loop-dependent part of the right-hand side
523            CALL WADD(N,Y,Z(1,istage),TMP)              ! TMP <- Y + Zi
524            CALL FUN_CHEM(T+rkC(istage)*H,TMP,RHS)      ! RHS <- Fun(Y+Zi)
525            ISTATUS(Nfun) = ISTATUS(Nfun) + 1
526!            RHS(1:N) = G(1:N) - Z(1:N,istage) + (H*rkGamma)*RHS(1:N)
527            CALL WSCAL(N, H*rkGamma, RHS, 1)
528            CALL WAXPY (N, -ONE, Z(1,istage), 1, RHS, 1)
529            CALL WAXPY (N, ONE, G,1, RHS,1)
530
531!~~~>   Solve the linear system
532            CALL SDIRK_Solve ( H, N, E, IP, IER, RHS )
533           
534!~~~>   Check convergence of Newton iterations
535            CALL SDIRK_ErrorNorm(N, RHS, SCAL, NewtonIncrement)
536            IF ( NewtonIter == 1 ) THEN
537                Theta      = ABS(ThetaMin)
538                NewtonRate = 2.0d0 
539            ELSE
540                Theta = NewtonIncrement/NewtonIncrementOld
541                IF (Theta < 0.99d0) THEN
542                    NewtonRate = Theta/(ONE-Theta)
543                    ! Predict error at the end of Newton process
544                    NewtonPredictedErr = NewtonIncrement &
545                               *Theta**(NewtonMaxit-NewtonIter)/(ONE-Theta)
546                    IF (NewtonPredictedErr >= NewtonTol) THEN
547                      ! Non-convergence of Newton: predicted error too large
548                      Qnewton = MIN(10.0d0,NewtonPredictedErr/NewtonTol)
549                      Fac = 0.8d0*Qnewton**(-ONE/(1+NewtonMaxit-NewtonIter))
550                      EXIT NewtonLoop
551                    END IF
552                ELSE ! Non-convergence of Newton: Theta too large
553                    EXIT NewtonLoop
554                END IF
555            END IF
556            NewtonIncrementOld = NewtonIncrement
557            ! Update solution: Z(:) <-- Z(:)+RHS(:)
558            CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1) 
559           
560            ! Check error in Newton iterations
561            NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol)
562            IF (NewtonDone) EXIT NewtonLoop
563           
564            END DO NewtonLoop
565           
566            IF (.NOT.NewtonDone) THEN
567                 !CALL RK_ErrorMsg(-12,T,H,Ierr);
568                 H = Fac*H; Reject=.TRUE.
569                 SkipJac = .TRUE.; SkipLU = .FALSE.
570                 CYCLE Tloop
571            END IF
572
573!~~~>  End of implified Newton iterations
574
575 
576   END DO stages
577
578
579!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
580!~~~>  Error estimation
581!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
582      ISTATUS(Nstp) = ISTATUS(Nstp) + 1
583      CALL Set2zero(N,TMP)
584      DO i = 1,rkS
585         IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,TMP,1)
586      END DO 
587
588      CALL SDIRK_Solve( H, N, E, IP, IER, TMP )
589      CALL SDIRK_ErrorNorm(N, TMP, SCAL, Err)
590
591!~~~> Computation of new step size Hnew
592      Fac  = FacSafe*(Err)**(-ONE/rkELO)
593      Fac  = MAX(FacMin,MIN(FacMax,Fac))
594      Hnew = H*Fac
595
596!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
597!~~~>  Accept/Reject step
598!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
599accept: IF ( Err < ONE ) THEN !~~~> Step is accepted
600
601         FirstStep=.FALSE.
602         ISTATUS(Nacc) = ISTATUS(Nacc) + 1
603
604!~~~> Update time and solution
605         T  =  T + H
606         ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:)
607         DO i = 1,rkS 
608            IF (rkD(i)/=ZERO) CALL WAXPY(N,rkD(i),Z(1,i),1,Y,1)
609         END DO 
610       
611!~~~> Update scaling coefficients
612         CALL SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL)
613
614!~~~> Next time step
615         Hnew = Tdirection*MIN(ABS(Hnew),Hmax)
616         ! Last T and H
617         RSTATUS(Ntexit) = T
618         RSTATUS(Nhexit) = H
619         RSTATUS(Nhnew)  = Hnew
620         ! No step increase after a rejection
621         IF (Reject) Hnew = Tdirection*MIN(ABS(Hnew),ABS(H))
622         Reject = .FALSE.
623         IF ((T+Hnew/Qmin-Tfinal)*Tdirection > ZERO) THEN
624            H = Tfinal-T
625         ELSE
626            Hratio=Hnew/H
627            ! If step not changed too much keep Jacobian and reuse LU
628            SkipLU = ( (Theta <= ThetaMin) .AND. (Hratio >= Qmin) &
629                                     .AND. (Hratio <= Qmax) )
630            IF (.NOT.SkipLU) H = Hnew
631         END IF
632         ! If convergence is fast enough, do not update Jacobian
633!         SkipJac = (Theta <= ThetaMin)
634         SkipJac = .FALSE.
635         
636      ELSE accept !~~~> Step is rejected
637
638         IF (FirstStep .OR. Reject) THEN
639             H = FacRej*H
640         ELSE
641             H = Hnew
642         END IF
643         Reject  = .TRUE.
644         SkipJac = .TRUE.
645         SkipLU  = .FALSE. 
646         IF (ISTATUS(Nacc) >= 1) ISTATUS(Nrej) = ISTATUS(Nrej) + 1
647         
648      END IF accept
649     
650      END DO Tloop
651
652      ! Successful return
653      Ierr  = 1
654 
655      END SUBROUTINE SDIRK_Integrator
656
657
658!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
659      SUBROUTINE SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL)
660!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
661      IMPLICIT NONE
662      INTEGER :: i, N, ITOL
663      KPP_REAL :: AbsTol(NVAR), RelTol(NVAR), &
664                       Y(NVAR), SCAL(NVAR)
665      IF (ITOL == 0) THEN
666        DO i=1,NVAR
667          SCAL(i) = ONE / ( AbsTol(1)+RelTol(1)*ABS(Y(i)) )
668        END DO
669      ELSE
670        DO i=1,NVAR
671          SCAL(i) = ONE / ( AbsTol(i)+RelTol(i)*ABS(Y(i)) )
672        END DO
673      END IF
674      END SUBROUTINE SDIRK_ErrorScale
675     
676     
677!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
678      SUBROUTINE SDIRK_ErrorNorm(N, Y, SCAL, Err)
679!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
680!     
681      INTEGER :: i, N
682      KPP_REAL :: Y(N), SCAL(N), Err     
683      Err = ZERO
684      DO i=1,N
685           Err = Err+(Y(i)*SCAL(i))**2
686      END DO
687      Err = MAX( SQRT(Err/DBLE(N)), 1.0d-10 )
688!
689      END SUBROUTINE SDIRK_ErrorNorm
690
691
692!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
693 SUBROUTINE SDIRK_ErrorMsg(Code,T,H,Ierr)
694!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
695!    Handles all error messages
696!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
697
698   KPP_REAL, INTENT(IN) :: T, H
699   INTEGER, INTENT(IN)  :: Code
700   INTEGER, INTENT(OUT) :: Ierr
701
702   Ierr = Code
703   PRINT * , &
704     'Forced exit from SDIRK due to the following error:'
705
706   SELECT CASE (Code)
707    CASE (-1)
708      PRINT * , '--> Improper value for maximal no of steps'
709    CASE (-2)
710      PRINT * , '--> Improper value for maximal no of Newton iterations'
711    CASE (-3)
712      PRINT * , '--> Hmin/Hmax/Hstart must be positive'
713    CASE (-4)
714      PRINT * , '--> FacMin/FacMax/FacRej must be positive'
715    CASE (-5)
716      PRINT * , '--> Improper tolerance values'
717    CASE (-6)
718      PRINT * , '--> No of steps exceeds maximum bound'
719    CASE (-7)
720      PRINT * , '--> Step size too small: T + 10*H = T', &
721            ' or H < Roundoff'
722    CASE (-8)
723      PRINT * , '--> Matrix is repeatedly singular'
724    CASE DEFAULT
725      PRINT *, 'Unknown Error code: ', Code
726   END SELECT
727
728   PRINT *, "T=", T, "and H=", H
729
730 END SUBROUTINE SDIRK_ErrorMsg
731     
732     
733!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
734      SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, &
735                   SkipJac, SkipLU, E, IP, Reject, ISING )
736!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
737!~~~>  Compute the matrix E = 1/(H*GAMMA)*Jac, and its decomposition
738!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
739     
740      IMPLICIT NONE
741     
742      KPP_REAL, INTENT(INOUT) :: H
743      KPP_REAL, INTENT(IN)    :: T, Y(NVAR)
744      LOGICAL, INTENT(INOUT)       :: SkipJac,SkipLU,Reject
745      INTEGER, INTENT(OUT)         :: ISING, IP(NVAR)
746#ifdef FULL_ALGEBRA
747      KPP_REAL, INTENT(INOUT) :: FJAC(NVAR,NVAR)
748      KPP_REAL, INTENT(OUT)   :: E(NVAR,NVAR)
749#else
750      KPP_REAL, INTENT(INOUT) :: FJAC(LU_NONZERO)
751      KPP_REAL, INTENT(OUT)   :: E(LU_NONZERO)
752#endif
753      KPP_REAL                :: HGammaInv
754      INTEGER                      :: i, j, ConsecutiveSng
755
756      ConsecutiveSng = 0
757      ISING = 1
758     
759Hloop: DO WHILE (ISING /= 0)
760     
761      HGammaInv = ONE/(H*rkGamma)
762
763!~~~>  Compute the Jacobian
764!      IF (SkipJac) THEN
765!          SkipJac = .FALSE.
766!      ELSE
767      IF (.NOT. SkipJac) THEN
768          CALL JAC_CHEM( T, Y, FJAC )
769          ISTATUS(Njac) = ISTATUS(Njac) + 1
770      END IF 
771     
772#ifdef FULL_ALGEBRA
773      DO j=1,NVAR
774         DO i=1,NVAR
775            E(i,j) = -FJAC(i,j)
776         END DO
777         E(j,j) = E(j,j)+HGammaInv
778      END DO
779      CALL DGETRF( NVAR, NVAR, E, NVAR, IP, ISING )
780#else
781      DO i = 1,LU_NONZERO
782            E(i) = -FJAC(i)
783      END DO
784      DO i = 1,NVAR
785          j = LU_DIAG(i); E(j) = E(j) + HGammaInv
786      END DO
787      CALL KppDecomp ( E, ISING)
788      IP(1) = 1
789#endif
790      ISTATUS(Ndec) = ISTATUS(Ndec) + 1
791
792      IF (ISING /= 0) THEN
793          WRITE (6,*) ' MATRIX IS SINGULAR, ISING=',ISING,' T=',T,' H=',H
794          ISTATUS(Nsng) = ISTATUS(Nsng) + 1; ConsecutiveSng = ConsecutiveSng + 1
795          IF (ConsecutiveSng >= 6) RETURN ! Failure
796          H = 0.5d0*H
797          SkipJac = .TRUE.
798          SkipLU  = .FALSE.
799          Reject  = .TRUE.
800      END IF
801     
802      END DO Hloop
803
804      END SUBROUTINE SDIRK_PrepareMatrix
805
806
807!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
808      SUBROUTINE SDIRK_Solve ( H, N, E, IP, ISING, RHS )
809!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
810!~~~>  Solves the system (H*Gamma-Jac)*x = RHS
811!      using the LU decomposition of E = I - 1/(H*Gamma)*Jac
812!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
813      IMPLICIT NONE
814      INTEGER, INTENT(IN)       :: N, IP(N), ISING
815      KPP_REAL, INTENT(IN) :: H
816#ifdef FULL_ALGEBRA
817      KPP_REAL, INTENT(IN) :: E(NVAR,NVAR)
818#else
819      KPP_REAL, INTENT(IN) :: E(LU_NONZERO)
820#endif
821      KPP_REAL, INTENT(INOUT) :: RHS(N)
822      KPP_REAL                :: HGammaInv
823     
824      HGammaInv = ONE/(H*rkGamma)
825      CALL WSCAL(N,HGammaInv,RHS,1)
826#ifdef FULL_ALGEBRA 
827      CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING )
828#else
829      CALL KppSolve(E, RHS)
830#endif
831      ISTATUS(Nsol) = ISTATUS(Nsol) + 1
832 
833      END SUBROUTINE SDIRK_Solve
834
835
836!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
837      SUBROUTINE Sdirk4a
838!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
839      sdMethod = S4A
840! Number of stages
841      rkS = 5
842
843! Method coefficients
844      rkGamma = .2666666666666666666666666666666667d0
845
846      rkA(1,1) = .2666666666666666666666666666666667d0
847      rkA(2,1) = .5000000000000000000000000000000000d0
848      rkA(2,2) = .2666666666666666666666666666666667d0
849      rkA(3,1) = .3541539528432732316227461858529820d0
850      rkA(3,2) = -.5415395284327323162274618585298197d-1
851      rkA(3,3) = .2666666666666666666666666666666667d0
852      rkA(4,1) = .8515494131138652076337791881433756d-1
853      rkA(4,2) = -.6484332287891555171683963466229754d-1
854      rkA(4,3) = .7915325296404206392428857585141242d-1
855      rkA(4,4) = .2666666666666666666666666666666667d0
856      rkA(5,1) = 2.100115700566932777970612055999074d0
857      rkA(5,2) = -.7677800284445976813343102185062276d0
858      rkA(5,3) = 2.399816361080026398094746205273880d0
859      rkA(5,4) = -2.998818699869028161397714709433394d0
860      rkA(5,5) = .2666666666666666666666666666666667d0
861
862      rkB(1)   = 2.100115700566932777970612055999074d0
863      rkB(2)   = -.7677800284445976813343102185062276d0
864      rkB(3)   = 2.399816361080026398094746205273880d0
865      rkB(4)   = -2.998818699869028161397714709433394d0
866      rkB(5)   = .2666666666666666666666666666666667d0
867
868      rkBhat(1)= 2.885264204387193942183851612883390d0
869      rkBhat(2)= -.1458793482962771337341223443218041d0
870      rkBhat(3)= 2.390008682465139866479830743628554d0
871      rkBhat(4)= -4.129393538556056674929560012190140d0
872      rkBhat(5)= 0.d0
873
874      rkC(1)   = .2666666666666666666666666666666667d0
875      rkC(2)   = .7666666666666666666666666666666667d0
876      rkC(3)   = .5666666666666666666666666666666667d0
877      rkC(4)   = .3661315380631796996374935266701191d0
878      rkC(5)   = 1.d0
879
880! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}
881      rkD(1)   = 0.d0
882      rkD(2)   = 0.d0
883      rkD(3)   = 0.d0
884      rkD(4)   = 0.d0
885      rkD(5)   = 1.d0
886
887! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}
888      rkE(1)   = -.6804000050475287124787034884002302d0
889      rkE(2)   = 1.558961944525217193393931795738823d0
890      rkE(3)   = -13.55893003128907927748632408763868d0
891      rkE(4)   = 15.48522576958521253098585004571302d0
892      rkE(5)   = 1.d0
893
894! Local order of Err estimate
895      rkElo    = 4
896
897! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}
898      rkTheta(2,1) = 1.875000000000000000000000000000000d0
899      rkTheta(3,1) = 1.708847304091539528432732316227462d0
900      rkTheta(3,2) = -.2030773231622746185852981969486824d0
901      rkTheta(4,1) = .2680325578937783958847157206823118d0
902      rkTheta(4,2) = -.1828840955527181631794050728644549d0
903      rkTheta(4,3) = .2968246986151577397160821594427966d0
904      rkTheta(5,1) = .9096171815241460655379433581446771d0
905      rkTheta(5,2) = -3.108254967778352416114774430509465d0
906      rkTheta(5,3) = 12.33727431701306195581826123274001d0
907      rkTheta(5,4) = -11.24557012450885560524143016037523d0
908
909! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
910      rkAlpha(2,1) = 2.875000000000000000000000000000000d0
911      rkAlpha(3,1) = .8500000000000000000000000000000000d0
912      rkAlpha(3,2) = .4434782608695652173913043478260870d0
913      rkAlpha(4,1) = .7352046091658870564637910527807370d0
914      rkAlpha(4,2) = -.9525565003057343527941920657462074d-1
915      rkAlpha(4,3) = .4290111305453813852259481840631738d0
916      rkAlpha(5,1) = -16.10898993405067684831655675112808d0
917      rkAlpha(5,2) = 6.559571569643355712998131800797873d0
918      rkAlpha(5,3) = -15.90772144271326504260996815012482d0
919      rkAlpha(5,4) = 25.34908987169226073668861694892683d0
920               
921!~~~> Coefficients for continuous solution
922!          rkD(1,1)= 24.74416644927758d0
923!          rkD(1,2)= -4.325375951824688d0
924!          rkD(1,3)= 41.39683763286316d0
925!          rkD(1,4)= -61.04144619901784d0
926!          rkD(1,5)= -3.391332232917013d0
927!          rkD(2,1)= -51.98245719616925d0
928!          rkD(2,2)= 10.52501981094525d0
929!          rkD(2,3)= -154.2067922191855d0
930!          rkD(2,4)= 214.3082125319825d0
931!          rkD(2,5)= 14.71166018088679d0
932!          rkD(3,1)= 33.14347947522142d0
933!          rkD(3,2)= -19.72986789558523d0
934!          rkD(3,3)= 230.4878502285804d0
935!          rkD(3,4)= -287.6629744338197d0
936!          rkD(3,5)= -18.99932366302254d0
937!          rkD(4,1)= -5.905188728329743d0
938!          rkD(4,2)= 13.53022403646467d0
939!          rkD(4,3)= -117.6778956422581d0
940!          rkD(4,4)= 134.3962081008550d0
941!          rkD(4,5)= 8.678995715052762d0
942!
943!         DO i=1,4  ! CONTi <-- Sum_j rkD(i,j)*Zj
944!           CALL Set2zero(N,CONT(1,i))
945!           DO j = 1,rkS
946!             CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1)
947!           END DO
948!         END DO
949         
950          rkELO = 4.0d0
951         
952      END SUBROUTINE Sdirk4a
953
954!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
955      SUBROUTINE Sdirk4b
956!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
957
958      sdMethod = S4B
959! Number of stages
960      rkS = 5
961
962! Method coefficients
963      rkGamma = .25d0
964
965      rkA(1,1) = 0.25d0
966      rkA(2,1) = 0.5d00
967      rkA(2,2) = 0.25d0
968      rkA(3,1) = 0.34d0
969      rkA(3,2) =-0.40d-1
970      rkA(3,3) = 0.25d0
971      rkA(4,1) = 0.2727941176470588235294117647058824d0
972      rkA(4,2) =-0.5036764705882352941176470588235294d-1
973      rkA(4,3) = 0.2757352941176470588235294117647059d-1
974      rkA(4,4) = 0.25d0
975      rkA(5,1) = 1.041666666666666666666666666666667d0
976      rkA(5,2) =-1.020833333333333333333333333333333d0
977      rkA(5,3) = 7.812500000000000000000000000000000d0
978      rkA(5,4) =-7.083333333333333333333333333333333d0
979      rkA(5,5) = 0.25d0
980
981      rkB(1)   =  1.041666666666666666666666666666667d0
982      rkB(2)   = -1.020833333333333333333333333333333d0
983      rkB(3)   =  7.812500000000000000000000000000000d0
984      rkB(4)   = -7.083333333333333333333333333333333d0
985      rkB(5)   =  0.250000000000000000000000000000000d0
986
987      rkBhat(1)=  1.069791666666666666666666666666667d0
988      rkBhat(2)= -0.894270833333333333333333333333333d0
989      rkBhat(3)=  7.695312500000000000000000000000000d0
990      rkBhat(4)= -7.083333333333333333333333333333333d0
991      rkBhat(5)=  0.212500000000000000000000000000000d0
992
993      rkC(1)   = 0.25d0
994      rkC(2)   = 0.75d0
995      rkC(3)   = 0.55d0
996      rkC(4)   = 0.50d0
997      rkC(5)   = 1.00d0
998
999! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}
1000      rkD(1)   = 0.0d0
1001      rkD(2)   = 0.0d0
1002      rkD(3)   = 0.0d0
1003      rkD(4)   = 0.0d0
1004      rkD(5)   = 1.0d0
1005
1006! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}
1007      rkE(1)   =  0.5750d0
1008      rkE(2)   =  0.2125d0
1009      rkE(3)   = -4.6875d0
1010      rkE(4)   =  4.2500d0
1011      rkE(5)   =  0.1500d0
1012
1013! Local order of Err estimate
1014      rkElo    = 4
1015
1016! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}
1017      rkTheta(2,1) = 2.d0
1018      rkTheta(3,1) = 1.680000000000000000000000000000000d0
1019      rkTheta(3,2) = -.1600000000000000000000000000000000d0
1020      rkTheta(4,1) = 1.308823529411764705882352941176471d0
1021      rkTheta(4,2) = -.1838235294117647058823529411764706d0
1022      rkTheta(4,3) = 0.1102941176470588235294117647058824d0
1023      rkTheta(5,1) = -3.083333333333333333333333333333333d0
1024      rkTheta(5,2) = -4.291666666666666666666666666666667d0
1025      rkTheta(5,3) =  34.37500000000000000000000000000000d0
1026      rkTheta(5,4) = -28.33333333333333333333333333333333d0
1027
1028! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1029      rkAlpha(2,1) = 3.
1030      rkAlpha(3,1) = .8800000000000000000000000000000000d0
1031      rkAlpha(3,2) = .4400000000000000000000000000000000d0
1032      rkAlpha(4,1) = .1666666666666666666666666666666667d0
1033      rkAlpha(4,2) = -.8333333333333333333333333333333333d-1
1034      rkAlpha(4,3) = .9469696969696969696969696969696970d0
1035      rkAlpha(5,1) = -6.d0
1036      rkAlpha(5,2) = 9.d0
1037      rkAlpha(5,3) = -56.81818181818181818181818181818182d0
1038      rkAlpha(5,4) = 54.d0
1039     
1040      END SUBROUTINE Sdirk4b
1041
1042!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1043      SUBROUTINE Sdirk2a
1044!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1045
1046      sdMethod = S2A
1047! Number of stages
1048      rkS = 2
1049
1050! Method coefficients
1051      rkGamma = .2928932188134524755991556378951510d0
1052
1053      rkA(1,1) = .2928932188134524755991556378951510d0
1054      rkA(2,1) = .7071067811865475244008443621048490d0
1055      rkA(2,2) = .2928932188134524755991556378951510d0
1056
1057      rkB(1)   = .7071067811865475244008443621048490d0
1058      rkB(2)   = .2928932188134524755991556378951510d0
1059
1060      rkBhat(1)= .6666666666666666666666666666666667d0
1061      rkBhat(2)= .3333333333333333333333333333333333d0
1062
1063      rkC(1)   = 0.292893218813452475599155637895151d0
1064      rkC(2)   = 1.0d0
1065
1066! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}
1067      rkD(1)   = 0.0d0
1068      rkD(2)   = 1.0d0
1069
1070! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}
1071      rkE(1)   =  0.4714045207910316829338962414032326d0
1072      rkE(2)   = -0.1380711874576983496005629080698993d0
1073
1074! Local order of Err estimate
1075      rkElo    = 2
1076
1077! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}
1078      rkTheta(2,1) = 2.414213562373095048801688724209698d0
1079
1080! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1081      rkAlpha(2,1) = 3.414213562373095048801688724209698d0
1082         
1083      END SUBROUTINE Sdirk2a
1084
1085!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1086      SUBROUTINE Sdirk2b
1087!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1088
1089      sdMethod = S2B
1090! Number of stages
1091      rkS      = 2
1092
1093! Method coefficients
1094      rkGamma  = 1.707106781186547524400844362104849d0
1095
1096      rkA(1,1) = 1.707106781186547524400844362104849d0
1097      rkA(2,1) = -.707106781186547524400844362104849d0
1098      rkA(2,2) = 1.707106781186547524400844362104849d0
1099
1100      rkB(1)   = -.707106781186547524400844362104849d0
1101      rkB(2)   = 1.707106781186547524400844362104849d0
1102
1103      rkBhat(1)= .6666666666666666666666666666666667d0
1104      rkBhat(2)= .3333333333333333333333333333333333d0
1105
1106      rkC(1)   = 1.707106781186547524400844362104849d0
1107      rkC(2)   = 1.0d0
1108
1109! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}
1110      rkD(1)   = 0.0d0
1111      rkD(2)   = 1.0d0
1112
1113! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}
1114      rkE(1)   = -.4714045207910316829338962414032326d0
1115      rkE(2)   =  .8047378541243650162672295747365659d0
1116
1117! Local order of Err estimate
1118      rkElo    = 2
1119
1120! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}
1121      rkTheta(2,1) = -.414213562373095048801688724209698d0
1122
1123! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1124      rkAlpha(2,1) = .5857864376269049511983112757903019d0
1125     
1126      END SUBROUTINE Sdirk2b
1127
1128
1129!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1130      SUBROUTINE Sdirk3a
1131!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1132
1133      sdMethod = S3A
1134! Number of stages
1135      rkS = 3
1136
1137! Method coefficients
1138      rkGamma = .2113248654051871177454256097490213d0
1139
1140      rkA(1,1) = .2113248654051871177454256097490213d0
1141      rkA(2,1) = .2113248654051871177454256097490213d0
1142      rkA(2,2) = .2113248654051871177454256097490213d0
1143      rkA(3,1) = .2113248654051871177454256097490213d0
1144      rkA(3,2) = .5773502691896257645091487805019573d0
1145      rkA(3,3) = .2113248654051871177454256097490213d0
1146
1147      rkB(1)   = .2113248654051871177454256097490213d0
1148      rkB(2)   = .5773502691896257645091487805019573d0
1149      rkB(3)   = .2113248654051871177454256097490213d0
1150
1151      rkBhat(1)= .2113248654051871177454256097490213d0
1152      rkBhat(2)= .6477918909913548037576239837516312d0
1153      rkBhat(3)= .1408832436034580784969504064993475d0
1154
1155      rkC(1)   = .2113248654051871177454256097490213d0
1156      rkC(2)   = .4226497308103742354908512194980427d0
1157      rkC(3)   = 1.d0
1158
1159! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}
1160      rkD(1)   = 0.d0
1161      rkD(2)   = 0.d0
1162      rkD(3)   = 1.d0
1163
1164! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}
1165      rkE(1)   =  0.9106836025229590978424821138352906d0
1166      rkE(2)   = -1.244016935856292431175815447168624d0
1167      rkE(3)   =  0.3333333333333333333333333333333333d0
1168
1169! Local order of Err estimate
1170      rkElo    = 2
1171
1172! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}
1173      rkTheta(2,1) =  1.0d0
1174      rkTheta(3,1) = -1.732050807568877293527446341505872d0
1175      rkTheta(3,2) =  2.732050807568877293527446341505872d0
1176
1177! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1178      rkAlpha(2,1) =   2.0d0
1179      rkAlpha(3,1) = -12.92820323027550917410978536602349d0
1180      rkAlpha(3,2) =   8.83012701892219323381861585376468d0
1181
1182      END SUBROUTINE Sdirk3a
1183
1184!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1185   END SUBROUTINE SDIRK ! AND ALL ITS INTERNAL PROCEDURES
1186!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1187
1188
1189!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1190      SUBROUTINE FUN_CHEM( T, Y, P )
1191!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1192
1193      USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO
1194      USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME
1195      USE KPP_ROOT_Function, ONLY: Fun
1196      USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1197      IMPLICIT NONE
1198
1199      KPP_REAL :: T, Told
1200      KPP_REAL :: Y(NVAR), P(NVAR)
1201     
1202      Told = TIME
1203      TIME = T
1204      CALL Update_SUN()
1205      CALL Update_RCONST()
1206     
1207      CALL Fun( Y, FIX, RCONST, P )
1208     
1209      TIME = Told
1210     
1211      END SUBROUTINE FUN_CHEM
1212
1213
1214!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1215      SUBROUTINE JAC_CHEM( T, Y, JV )
1216!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1217
1218      USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO
1219      USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME
1220      USE KPP_ROOT_Jacobian
1221      USE KPP_ROOT_Jacobian, ONLY: Jac_SP
1222      USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1223      IMPLICIT NONE
1224 
1225      KPP_REAL ::  T, Told
1226      KPP_REAL ::  Y(NVAR)
1227#ifdef FULL_ALGEBRA
1228      KPP_REAL :: JS(LU_NONZERO), JV(NVAR,NVAR)
1229      INTEGER :: i, j
1230#else
1231      KPP_REAL :: JV(LU_NONZERO)
1232#endif
1233 
1234      Told = TIME
1235      TIME = T
1236      CALL Update_SUN()
1237      CALL Update_RCONST()
1238
1239#ifdef FULL_ALGEBRA
1240      CALL Jac_SP(Y, FIX, RCONST, JS)
1241      DO j=1,NVAR
1242         DO i=1,NVAR
1243          JV(i,j) = 0.0D0
1244         END DO
1245      END DO
1246      DO i=1,LU_NONZERO
1247         JV(LU_IROW(i),LU_ICOL(i)) = JS(i)
1248      END DO
1249#else
1250      CALL Jac_SP(Y, FIX, RCONST, JV)
1251#endif
1252      TIME = Told
1253
1254      END SUBROUTINE JAC_CHEM
1255
1256END MODULE KPP_ROOT_Integrator
1257
1258
Note: See TracBrowser for help on using the repository browser.