source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/rosenbrock_vec.f90 @ 2696

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

Merge of branch palm4u into trunk

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