C----------------------------------------------------------------------- C Demonstration program for SVODPK. C ODE system from ns-species interaction PDE in 2 dimensions. C C This is the version of 23 April 2002. C C This version is in single precision. C----------------------------------------------------------------------- C This program solves a stiff ODE system that arises from a system C of partial differential equations. The PDE system is a food web C population model, with predator-prey interaction and diffusion on C the unit square in two dimensions. The dependent variable vector is C 1 2 ns C c = (c , c , ..., c ) C and the PDEs are as follows: C i i i C dc /dt = d(i)*(c + c ) + f (x,y,c) (i=1,...,ns) C xx yy i C where C i ns j C f (x,y,c) = c *(b(i) + sum a(i,j)*c ) . C i j=1 C The number of species is ns = 2*np, with the first np being prey and C the last np being predators. The coefficients a(i,j), b(i), d(i) are C a(i,i) = -a (all i) C a(i,j) = -g (i .le. np, j .gt. np) C a(i,j) = e (i .gt. np, j .le. np) C b(i) = b*(1 + alpha*x*y) (i .le. np) C b(i) = -b*(1 + alpha*x*y) (i .gt. np) C d(i) = dprey (i .le. np) C d(i) = dpred (i .gt. np) C The various scalar parameters are set in Subroutine SETPAR. C C The boundary conditions are: normal derivative = 0. C A polynomial in x and y is used to set the initial conditions. C C The PDEs are discretized by central differencing on a mx by my mesh. C C The ODE system is solved by SVODPK using method flag values C MF = 10, 21, and 29. The final time is TMAX = 10, except C that for MF = 10 it is TMAX = 1.0e-3 because the problem is stiff. C C Two preconditioner matrices are used. One uses a fixed number of C Gauss-Seidel iterations based on the diffusion terms only. C The other preconditioner is a block-diagonal matrix based on C the partial derivatives of the interaction terms f only, using C block-grouping (computing only a subset of the ns by ns blocks). C For MF = 21, these two preconditioners are applied on C the left and right, respectively. C For MF = 29, the inverse of the product is applied. C C Two output files are written: one with the problem description and C and performance statistics on unit 6, and one with solution profiles C at selected output times (for MF = 21 only) on unit 8. C----------------------------------------------------------------------- C Note: In addition to the main program and 10 subroutines given below, C this program requires the LINPACK subroutines SGEFA and SGESL, and C the BLAS routine SAXPY. C----------------------------------------------------------------------- C Reference: C Peter N. Brown and Alan C. Hindmarsh, C Reduced Storage Matrix Methods in Stiff ODE Systems, C J. Appl. Math. & Comp., 31 (1989), pp. 40-91. C Also available as Lawrence Livermore National Laboratory C Report UCRL-95088, Rev. 1, June 1987. C----------------------------------------------------------------------- EXTERNAL FWEB, JACBG, SOLSBG INTEGER NS, MX, MY, MXNS, 1 MP, MQ, MPSQ, ITMAX, 2 MESHX,MESHY,NGX,NGY,NGRP,MXMP,JGX,JGY,JIGX,JIGY,JXR,JYR INTEGER I, IMOD3, IOPT, IPAR, IOUT, ISTATE, ITASK, ITOL, IWORK, 1 JACFLG, JPRE, LENIW, LENRW, LIW, LRW, MF, NCFL, NCFN, 2 NETF, NEQ, NFE, NFLDIF, NFNDIF, NLI, NLIDIF, NNI, NNIDIF, 3 NOUT, NPE, NPS, NQU, NSDIF, NST REAL AA, EE, GG, BB, DPREY, DPRED, 1 AX, AY, ACOEF, BCOEF, DX, DY, ALPH, DIFF, COX, COY, 2 UROUND, SRUR REAL AVDIM, ATOL, CC, HU, RCFL, RCFN, RUMACH, 1 RPAR, RTOL, RWORK, T, TOUT C C The problem COMMON blocks below allow for up to 20 species, C up to a 50x50 mesh, and up to a 20x20 group structure. COMMON /PCOM0/ AA, EE, GG, BB, DPREY, DPRED COMMON /PCOM1/ NS, MX, MY, MXNS, AX, AY, ACOEF(20,20), BCOEF(20), 1 DX, DY, ALPH, DIFF(20), COX(20), COY(20) COMMON /PCOM2/ MP, MQ, MPSQ, ITMAX, UROUND, SRUR COMMON /PCOM3/ MESHX, MESHY, NGX, NGY, NGRP, MXMP, 2 JGX(21), JGY(21), JIGX(50), JIGY(50), JXR(20), JYR(20) C C The dimension of CC and RPAR below must be .ge. NEQ = NS*MX*MY. C The dimension LRW of RWORK must be .ge. 17*NEQ + NS*NS*NGRP + 61, C and the dimension LIW of IWORK must be .ge. 30 + NS*NGRP. DIMENSION CC(216), RWORK(3877), IWORK(54), RPAR(216) DATA LRW/3877/, LIW/54/ C OPEN (unit=6, file='demout', status='new') OPEN (unit=8, file='ccout', status='new') C AX = 1.0E0 AY = 1.0E0 C C Call setpar to set problem parameters. CALL SETPAR C WRITE(6,10)NS 10 FORMAT(' Demonstration program for SVODPK package'// 1 ' Food web problem with NS species, NS =',I4/ 2 ' Predator-prey interaction and diffusion on a 2-D square'//) WRITE(6,14) AA,EE,GG,BB,DPREY,DPRED 14 FORMAT(' Matrix parameters: a =',E12.4,' e =',E12.4, 1 ' g =',E12.4//' b parameter =',E12.4// 2 ' Diffusion coefficients: Dprey =',E12.4,' Dpred =',E12.4/) WRITE(6,16) ALPH 16 FORMAT(' Rate parameter alpha =',E12.4//) C C Set remaining problem parameters. NEQ = NS*MX*MY MXNS = MX*NS DX = AX/(MX-1) DY = AY/(MY-1) DO 20 I = 1,NS COX(I) = DIFF(I)/DX**2 COY(I) = DIFF(I)/DY**2 20 CONTINUE C WRITE(6,25) MX,MY,NEQ 25 FORMAT(' Mesh dimensions (mx,my) =',2I4, 1 //' Total system size is NEQ =',I7///) C C Set remaining method parameters. JPRE = 3 JACFLG = 1 IWORK(3) = JPRE IWORK(4) = JACFLG IOPT = 0 MP = NS MQ = MX*MY MPSQ = NS*NS UROUND = RUMACH() SRUR = SQRT(UROUND) MESHX = MX MESHY = MY MXMP = MESHX*MP NGX = 2 NGY = 2 NGRP = NGX*NGY CALL GSET (MESHX, NGX, JGX, JIGX, JXR) CALL GSET (MESHY, NGY, JGY, JIGY, JYR) IWORK(1) = MPSQ*NGRP IWORK(2) = MP*NGRP ITMAX = 5 ITOL = 1 RTOL = 1.0E-4 ATOL = RTOL ITASK = 1 WRITE(6,30)NGRP,NGX,NGY,ITMAX,RTOL,ATOL 30 FORMAT(' Preconditioning uses interaction-only block-diagonal', 1 ' matrix'/' with block-grouping, and G-S iterations'// 2 ' Number of diagonal block groups = ngrp =',I4, 3 ' (ngx by ngy, ngx =',I2,' ngy =',I2,' )'// 4 ' G-S preconditioner uses itmax iterations, itmax =',I3/// 5 ' Tolerance parameters: RTOL =',E10.2,' ATOL =',E10.2//) C C C Loop over MF values 10, 21, 29. C DO 90 MF = 10,29 IF (MF .GT. 10 .AND. MF .LT. 21) GO TO 90 IF (MF .GT. 21 .AND. MF .LT. 29) GO TO 90 WRITE(6,40)MF 40 FORMAT(85('*')//' Solution with MF =',I3// 1 ' t',7X,'NSTEP',' NFE',' NNI',' NLI',' NPE',' NQ', 2 4X,'H',10X,'AVDIM',' NCF rate LCF rate') C T = 0.0E0 TOUT = 1.0E-8 NOUT = 18 IF (MF .EQ. 10) NOUT = 6 CALL CINIT (CC) IF (MF .EQ. 21) CALL OUTWEB (T, CC, NS, MX, MY, 8) ISTATE = 1 NLI = 0 NNI = 0 NCFN = 0 NCFL = 0 NST = 0 C C Loop over output times and call SVODPK. C DO 70 IOUT = 1,NOUT CALL SVODPK (FWEB, NEQ, CC, T, TOUT, ITOL, RTOL, ATOL, ITASK, 1 ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JACBG, SOLSBG, 2 MF, RPAR, IPAR) NSDIF = IWORK(11) - NST NST = IWORK(11) NFE = IWORK(12) NPE = IWORK(13) NNIDIF = IWORK(20) - NNI NNI = IWORK(20) NLIDIF = IWORK(23) - NLI NLI = IWORK(23) NFNDIF = IWORK(21) - NCFN NCFN = IWORK(21) NFLDIF = IWORK(25) - NCFL NCFL = IWORK(25) NQU = IWORK(14) HU = RWORK(11) AVDIM = 0.0E0 RCFN = 0.0E0 RCFL = 0.0E0 IF (NNIDIF .GT. 0) AVDIM = REAL(NLIDIF)/REAL(NNIDIF) IF (NSDIF .GT. 0) RCFN = REAL(NFNDIF)/REAL(NSDIF) IF (NNIDIF .GT. 0) RCFL = REAL(NFLDIF)/REAL(NNIDIF) WRITE(6,50)T,NST,NFE,NNI,NLI,NPE,NQU,HU,AVDIM,RCFN,RCFL 50 FORMAT(E10.2,I5,I6,3I5,I4,2E11.2,E10.2,E12.2) IMOD3 = IOUT - 3*(IOUT/3) IF (MF .EQ. 21 .AND. IMOD3 .EQ. 0) CALL OUTWEB (T,CC,NS,MX,MY,8) IF (ISTATE .EQ. 2) GO TO 65 WRITE(6,60)T 60 FORMAT(//' Final time reached =',E12.4//) GO TO 75 65 CONTINUE IF (TOUT .GT. 0.9E0) TOUT = TOUT + 1.0E0 IF (TOUT .LT. 0.9E0) TOUT = TOUT*10.0E0 70 CONTINUE C 75 CONTINUE NST = IWORK(11) NFE = IWORK(12) NPE = IWORK(13) LENRW = IWORK(17) LENIW = IWORK(18) NNI = IWORK(20) NLI = IWORK(23) NPS = IWORK(24) IF (NNI .GT. 0) AVDIM = REAL(NLI)/REAL(NNI) NCFN = IWORK(21) NCFL = IWORK(25) NETF = IWORK(22) WRITE (6,80) LENRW,LENIW,NST,NFE,NPE,NPS,NNI,NLI,AVDIM, 1 NCFN,NCFL,NETF 80 FORMAT(//' Final statistics for this run:'/ 1 ' RWORK size =',I8,' IWORK size =',I6/ 2 ' Number of time steps =',I5/ 3 ' Number of f evaluations =',I5/ 4 ' Number of preconditioner evals. =',I5/ 4 ' Number of preconditioner solves =',I5/ 5 ' Number of nonlinear iterations =',I5/ 5 ' Number of linear iterations =',I5/ 6 ' Average subspace dimension =',F8.4/ 7 I5,' nonlinear conv. failures,',I5,' linear conv. failures,', 8 I5,' error test failures'//) C 90 CONTINUE STOP C------ End of Main Program for SVODPK demonstration program ----------- END SUBROUTINE SETPAR C----------------------------------------------------------------------- C This routine sets the problem parameters. C It sets NS, MX, MY, and problem coefficients ACOEF, BCOEF, DIFF, C and ALPH, using parameters NP, AA, EE, GG, BB, DPREY, DPRED. C----------------------------------------------------------------------- INTEGER NS, MX, MY, MXNS INTEGER I, J, NP REAL AA, EE, GG, BB, DPREY, DPRED, 1 AX, AY, ACOEF, BCOEF, DX, DY, ALPH, DIFF, COX, COY COMMON /PCOM0/ AA, EE, GG, BB, DPREY, DPRED COMMON /PCOM1/ NS, MX, MY, MXNS, AX, AY, ACOEF(20,20), BCOEF(20), 1 DX, DY, ALPH, DIFF(20), COX(20), COY(20) C NP = 3 MX = 6 MY = 6 AA = 1.0E0 EE = 1.0E4 GG = 0.5E-6 BB = 1.0E0 DPREY = 1.0E0 DPRED = 0.5E0 ALPH = 1.0E0 NS = 2*NP DO 70 J = 1,NP DO 60 I = 1,NP ACOEF(NP+I,J) = EE ACOEF(I,NP+J) = -GG 60 CONTINUE ACOEF(J,J) = -AA ACOEF(NP+J,NP+J) = -AA BCOEF(J) = BB BCOEF(NP+J) = -BB DIFF(J) = DPREY DIFF(NP+J) = DPRED 70 CONTINUE C RETURN C------------ End of Subroutine SETPAR ------------------------------- END SUBROUTINE GSET (M, NG, JG, JIG, JR) C----------------------------------------------------------------------- C This routine sets arrays JG, JIG, and JR describing C a uniform partition of (1,2,...,M) into NG groups. C----------------------------------------------------------------------- INTEGER M, NG, JG, JIG, JR DIMENSION JG(*), JIG(*), JR(*) INTEGER IG, J, LEN1, MPER, NGM1 C MPER = M/NG DO 10 IG = 1,NG 10 JG(IG) = 1 + (IG - 1)*MPER JG(NG+1) = M + 1 C NGM1 = NG - 1 LEN1 = NGM1*MPER DO 20 J = 1,LEN1 20 JIG(J) = 1 + (J-1)/MPER LEN1 = LEN1 + 1 DO 25 J = LEN1,M 25 JIG(J) = NG C DO 30 IG = 1,NGM1 30 JR(IG) = 0.5E0 + (IG - 0.5E0)*MPER JR(NG) = 0.5E0*(1 + NGM1*MPER + M) C RETURN C------------ End of Subroutine GSET --------------------------------- END SUBROUTINE CINIT (CC) C----------------------------------------------------------------------- C This routine computes and loads the vector of initial values. C----------------------------------------------------------------------- REAL CC DIMENSION CC(*) INTEGER NS, MX, MY, MXNS INTEGER I, ICI, IOFF, IYOFF, JX, JY REAL AX,AY,ACOEF,BCOEF,DX,DY,ALPH,DIFF,COX,COY REAL ARGX, ARGY, X, Y COMMON /PCOM1/ NS, MX, MY, MXNS, AX, AY, ACOEF(20,20), BCOEF(20), 1 DX, DY, ALPH, DIFF(20), COX(20), COY(20) C DO 20 JY = 1,MY Y = (JY-1)*DY ARGY = 16.0E0*Y*Y*(AY-Y)*(AY-Y) IYOFF = MXNS*(JY-1) DO 10 JX = 1,MX X = (JX-1)*DX ARGX = 16.0E0*X*X*(AX-X)*(AX-X) IOFF = IYOFF + NS*(JX-1) DO 5 I = 1,NS ICI = IOFF + I CC(ICI) = 10.0E0 + I*ARGX*ARGY 5 CONTINUE 10 CONTINUE 20 CONTINUE RETURN C------------ End of Subroutine CINIT -------------------------------- END SUBROUTINE OUTWEB (T, C, NS, MX, MY, LUN) C----------------------------------------------------------------------- C This routine prints the values of the individual species densities C at the current time T. All output is on unit LUN. C----------------------------------------------------------------------- INTEGER NS, MX, MY, LUN REAL T, C DIMENSION C(NS,MX,MY) INTEGER I, JX, JY C WRITE(LUN,10) T 10 FORMAT(/80('-')/30X,'At time t = ',E16.8/80('-') ) C DO 40 I = 1,NS WRITE(LUN,20) I 20 FORMAT(' The species c(',I2,') values are') DO 30 JY = MY,1,-1 WRITE(LUN,25) (C(I,JX,JY),JX=1,MX) 25 FORMAT(6(G13.6)) 30 CONTINUE WRITE(LUN,35) 35 FORMAT(80('-'),/) 40 CONTINUE C RETURN C------------ End of Subroutine OUTWEB ------------------------------- END SUBROUTINE FWEB (NEQ, T, CC, CDOT, RPAR, IPAR) C----------------------------------------------------------------------- C This routine computes the right-hand side of the ODE system and C returns it in CDOT. C The interaction rates are computed by calls to WEBR, and these are C saved in RPAR for use in preconditioning. C----------------------------------------------------------------------- INTEGER NEQ, IPAR REAL T, CC, CDOT, RPAR DIMENSION CC(*), CDOT(*), RPAR(*) INTEGER NS, MX, MY, MXNS INTEGER I, IC, ICI, IDXL, IDXU, IDYL, IDYU, IYOFF, JX, JY REAL AX,AY,ACOEF,BCOEF,DX,DY,ALPH,DIFF,COX,COY REAL DCXLI, DCXUI, DCYLI, DCYUI, X, Y COMMON /PCOM1/ NS, MX, MY, MXNS, AX, AY, ACOEF(20,20), BCOEF(20), 1 DX, DY, ALPH, DIFF(20), COX(20), COY(20) C DO 100 JY = 1,MY Y = (JY-1)*DY IYOFF = MXNS*(JY-1) IDYU = MXNS IF (JY .EQ. MY) IDYU = -MXNS IDYL = MXNS IF (JY .EQ. 1) IDYL = -MXNS DO 90 JX = 1,MX X = (JX-1)*DX IC = IYOFF + NS*(JX-1) + 1 C Get interaction rates at one point (x,y). CALL WEBR (X, Y, T, CC(IC), RPAR(IC)) IDXU = NS IF (JX .EQ. MX) IDXU = -NS IDXL = NS IF (JX .EQ. 1) IDXL = -NS DO 80 I = 1,NS ICI = IC + I - 1 C Do differencing in y. DCYLI = CC(ICI) - CC(ICI-IDYL) DCYUI = CC(ICI+IDYU) - CC(ICI) C Do differencing in x. DCXLI = CC(ICI) - CC(ICI-IDXL) DCXUI = CC(ICI+IDXU) - CC(ICI) C Collect terms and load CDOT elements. CDOT(ICI) = COY(I)*(DCYUI - DCYLI) + COX(I)*(DCXUI - DCXLI) 1 + RPAR(ICI) 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN C------------ End of Subroutine FWEB --------------------------------- END SUBROUTINE WEBR (X, Y, T, C, RATE) C----------------------------------------------------------------------- C This routine computes the interaction rates for the species C C(1),...,C(NS), at one spatial point and at time T. C----------------------------------------------------------------------- REAL X, Y, T, C, RATE DIMENSION C(*), RATE(*) INTEGER NS, MX, MY, MXNS INTEGER I REAL AX,AY,ACOEF,BCOEF,DX,DY,ALPH,DIFF,COX,COY REAL FAC COMMON /PCOM1/ NS, MX, MY, MXNS, AX, AY, ACOEF(20,20), BCOEF(20), 1 DX, DY, ALPH, DIFF(20), COX(20), COY(20) C DO 10 I = 1,NS 10 RATE(I) = 0.0E0 DO 15 I = 1,NS CALL SAXPY (NS, C(I), ACOEF(1,I), 1, RATE, 1) 15 CONTINUE FAC = 1.0E0 + ALPH*X*Y DO 20 I = 1,NS 20 RATE(I) = C(I)*(BCOEF(I)*FAC + RATE(I)) RETURN C------------ End of Subroutine WEBR --------------------------------- END SUBROUTINE JACBG (F, NEQ, T, CC, CCSV, REWT, F0, F1, HB0, 1 BD, IPBD, IER, RPAR, IPAR) C----------------------------------------------------------------------- C This routine generates part of the block-diagonal part of the C Jacobian, multiplies by -HB0, adds the identity matrix, C and calls SGEFA to do LU decomposition of each diagonal block. C The computation of the diagonal blocks uses the block and grouping C information in /PCOM1/ and /PCOM2/. One block per group is computed. C The Jacobian elements are generated by difference quotients C using calls to the routine FBG. C----------------------------------------------------------------------- C The two COMMON blocks below are used for internal communication. C The variables used are: C MP = size of blocks in block-diagonal preconditioning matrix. C MQ = number of blocks in each direction (NEQ = MP*MQ). C MPSQ = MP*MP. C UROUND = unit roundoff, generated by a call UROUND = RUMACH(). C SRUR = SQRT(UROUND). C MESHX = x mesh size C MESHY = y mesh size (mesh is MESHX by MESHY) C NGX = no. groups in x direction in block-grouping scheme. C NGY = no. groups in y direction in block-grouping scheme. C NGRP = total number of groups = NGX*NGY. C MXMP = MESHX*MP. C JGX = length NGX+1 array of group boundaries in x direction. C group IGX has x indices JX = JGX(IGX),...,JGX(IGX+1)-1. C JIGX = length MESHX array of x group indices vs x node index. C x node index JX is in x group JIGX(JX). C JXR = length NGX array of x indices representing the x groups. C the index for x group IGX is JX = JXR(IGX). C JGY, JIGY, JYR = analogous arrays for grouping in y direction. C----------------------------------------------------------------------- EXTERNAL F INTEGER NEQ, IPBD, IER, IPAR REAL T, CC, CCSV, REWT, F0, F1, HB0, BD, RPAR DIMENSION CC(*), CCSV(*), REWT(*), F0(*), F1(*), BD(*), IPBD(*), 1 RPAR(*) INTEGER MP, MQ, MPSQ, ITMAX, 2 MESHX,MESHY,NGX,NGY,NGRP,MXMP,JGX,JGY,JIGX,JIGY,JXR,JYR INTEGER I, IBD, IDIAG, IF0, IF00, IG, IGX, IGY, IIP, 1 J, JJ, JX, JY REAL UROUND, SRUR REAL FAC, R, R0, SVNORM C COMMON /PCOM2/ MP, MQ, MPSQ, ITMAX, UROUND, SRUR COMMON /PCOM3/ MESHX, MESHY, NGX, NGY, NGRP, MXMP, 1 JGX(21), JGY(21), JIGX(50), JIGY(50), JXR(20), JYR(20) C C----------------------------------------------------------------------- C Make MP calls to FBG to approximate each diagonal block of Jacobian. C Here, RPAR contains the base fb value. C R0 is a minimum increment factor for the difference quotient. C----------------------------------------------------------------------- FAC = SVNORM (NEQ, F0, REWT) R0 = 1000.0E0*ABS(HB0)*UROUND*NEQ*FAC IF (R0 .EQ. 0.0E0) R0 = 1.0E0 IBD = 0 DO 240 IGY = 1,NGY JY = JYR(IGY) IF00 = (JY - 1)*MXMP DO 230 IGX = 1,NGX JX = JXR(IGX) IF0 = IF00 + (JX - 1)*MP DO 220 J = 1,MP JJ = IF0 + J R = MAX(SRUR*ABS(CC(JJ)),R0/REWT(JJ)) CC(JJ) = CC(JJ) + R FAC = -HB0/R CALL FBG (T, CC, JX, JY, F1) DO 210 I = 1,MP 210 BD(IBD+I) = (F1(I) - RPAR(IF0+I))*FAC CC(JJ) = CCSV(JJ) IBD = IBD + MP 220 CONTINUE 230 CONTINUE 240 CONTINUE C C Add identity matrix and do LU decompositions on blocks. 260 CONTINUE IBD = 1 IIP = 1 DO 280 IG = 1,NGRP IDIAG = IBD DO 270 I = 1,MP BD(IDIAG) = BD(IDIAG) + 1.0E0 270 IDIAG = IDIAG + (MP + 1) CALL SGEFA (BD(IBD), MP, MP, IPBD(IIP), IER) IF (IER .NE. 0) GO TO 290 IBD = IBD + MPSQ IIP = IIP + MP 280 CONTINUE 290 RETURN C------------ End of Subroutine JACBG -------------------------------- END SUBROUTINE FBG (T, CC, JX, JY, CDOT) C----------------------------------------------------------------------- C This routine computes one block of the interaction terms of the C system, namely block (JX,JY), for use in preconditioning. C----------------------------------------------------------------------- INTEGER JX, JY REAL T, CC, CDOT DIMENSION CC(*), CDOT(*) INTEGER NS, MX, MY, MXNS INTEGER IBLOK, IC REAL AX,AY,ACOEF,BCOEF,DX,DY,ALPH,DIFF,COX,COY REAL X, Y C COMMON /PCOM1/ NS, MX, MY, MXNS, AX, AY, ACOEF(20,20), BCOEF(20), 1 DX, DY, ALPH, DIFF(20), COX(20), COY(20) C IBLOK = JX + (JY-1)*MX Y = (JY-1)*DY X = (JX-1)*DX IC = NS*(IBLOK-1) + 1 CALL WEBR (X, Y, T, CC(IC), CDOT) RETURN C------------ End of Subroutine FBG ---------------------------------- END SUBROUTINE SOLSBG (N, T, CC, F0, WK, HB0, BD, IPBD, V, LR, 1 IER, RPAR, IPAR) C----------------------------------------------------------------------- C This routine applies one of two inverse preconditioner matrices C to the array V, using the interaction-only block-diagonal Jacobian C with block-grouping, denoted Jr, and Gauss-Seidel applied to the C diffusion contribution to the Jacobian, denoted Jd. C When LR = 1, it calls GS for a Gauss-Seidel approximation C to ((I-HB0*Jd)-inverse)*V, and stores the result in V. C when LR = 2, it computes ((I-HB0*Jr)-inverse)*V, using LU C factors of the blocks in BD, and pivot information in IPBD. C In both cases, the array V is overwritten with the solution vector. C----------------------------------------------------------------------- INTEGER N, IPBD, LR, IER, IPAR REAL T, CC, F0, WK, HB0, BD, V, RPAR DIMENSION CC(*), F0(*), WK(*), BD(*), IPBD(*), V(*) INTEGER MP, MQ, MPSQ, ITMAX, 2 MESHX,MESHY,NGX,NGY,NGRP,MXMP,JGX,JGY,JIGX,JIGY,JXR,JYR INTEGER IBD, IG0, IGM1, IGX, IGY, IIP, IV, JX, JY REAL UROUND, SRUR C COMMON /PCOM2/ MP, MQ, MPSQ, ITMAX, UROUND, SRUR COMMON /PCOM3/ MESHX, MESHY, NGX, NGY, NGRP, MXMP, 1 JGX(21), JGY(21), JIGX(50), JIGY(50), JXR(20), JYR(20) C IER = 0 C IF (LR .EQ. 1) THEN CALL GS (N, HB0, V, WK) ENDIF IF (LR .EQ. 2) THEN IV = 1 DO 20 JY = 1,MESHY IGY = JIGY(JY) IG0 = (IGY - 1)*NGX DO 10 JX = 1,MESHX IGX = JIGX(JX) IGM1 = IGX - 1 + IG0 IBD = 1 + IGM1*MPSQ IIP = 1 + IGM1*MP CALL SGESL (BD(IBD), MP, MP, IPBD(IIP), V(IV), 0) IV = IV + MP 10 CONTINUE 20 CONTINUE ENDIF C RETURN C------------ End of Subroutine SOLSBG ------------------------------- END SUBROUTINE GS (N, HB0, Z, X) C----------------------------------------------------------------------- C This routine performs ITMAX Gauss-Seidel iterations C to compute an approximation to P-inverse*Z, where P = I - HB0*Jd, and C Jd represents the diffusion contributions to the Jacobian. C Z contains the answer on return. C The dimensions below assume NS .le. 20. C----------------------------------------------------------------------- INTEGER N REAL HB0, Z, X DIMENSION Z(*), X(*) INTEGER NS, MX, MY, MXNS, 1 MP, MQ, MPSQ, ITMAX INTEGER I, IC, ICI, ITER, IYOFF, JX, JY REAL AX,AY,ACOEF,BCOEF,DX,DY,ALPH,DIFF,COX,COY, 2 UROUND, SRUR REAL BETA,BETA2,COF1,ELAMDA,GAMMA,GAMMA2 DIMENSION BETA(20), GAMMA(20), BETA2(20), GAMMA2(20), COF1(20) COMMON /PCOM1/ NS, MX, MY, MXNS, AX, AY, ACOEF(20,20), BCOEF(20), 1 DX, DY, ALPH, DIFF(20), COX(20), COY(20) COMMON /PCOM2/ MP, MQ, MPSQ, ITMAX, UROUND, SRUR C C----------------------------------------------------------------------- C Write matrix as P = D - L - U. C Load local arrays BETA, BETA2, GAMMA, GAMMA2, and COF1. C----------------------------------------------------------------------- DO 10 I = 1,NS ELAMDA = 1.E0/(1.E0 + 2.E0*HB0*(COX(I) + COY(I))) BETA(I) = HB0*COX(I)*ELAMDA BETA2(I) = 2.E0*BETA(I) GAMMA(I) = HB0*COY(I)*ELAMDA GAMMA2(I) = 2.E0*GAMMA(I) COF1(I) = ELAMDA 10 CONTINUE C----------------------------------------------------------------------- C Begin iteration loop. C Load array X with (D-inverse)*Z for first iteration. C----------------------------------------------------------------------- ITER = 1 C DO 50 JY = 1,MY IYOFF = MXNS*(JY-1) DO 40 JX = 1,MX IC = IYOFF + NS*(JX-1) DO 30 I = 1,NS ICI = IC + I X(ICI) = COF1(I)*Z(ICI) Z(ICI) = 0.E0 30 CONTINUE 40 CONTINUE 50 CONTINUE GO TO 160 C----------------------------------------------------------------------- C Calculate (D-inverse)*U*X. C----------------------------------------------------------------------- 70 CONTINUE ITER = ITER + 1 JY = 1 JX = 1 IC = NS*(JX-1) DO 75 I = 1,NS ICI = IC + I 75 X(ICI) = BETA2(I)*X(ICI+NS) + GAMMA2(I)*X(ICI+MXNS) DO 85 JX = 2,MX-1 IC = NS*(JX-1) DO 80 I = 1,NS ICI = IC + I 80 X(ICI) = BETA(I)*X(ICI+NS) + GAMMA2(I)*X(ICI+MXNS) 85 CONTINUE JX = MX IC = NS*(JX-1) DO 90 I = 1,NS ICI = IC + I 90 X(ICI) = GAMMA2(I)*X(ICI+MXNS) DO 115 JY = 2,MY-1 IYOFF = MXNS*(JY-1) JX = 1 IC = IYOFF DO 95 I = 1,NS ICI = IC + I 95 X(ICI) = BETA2(I)*X(ICI+NS) + GAMMA(I)*X(ICI+MXNS) DO 105 JX = 2,MX-1 IC = IYOFF + NS*(JX-1) DO 100 I = 1,NS ICI = IC + I 100 X(ICI) = BETA(I)*X(ICI+NS) + GAMMA(I)*X(ICI+MXNS) 105 CONTINUE JX = MX IC = IYOFF + NS*(JX-1) DO 110 I = 1,NS ICI = IC + I 110 X(ICI) = GAMMA(I)*X(ICI+MXNS) 115 CONTINUE JY = MY IYOFF = MXNS*(JY-1) JX = 1 IC = IYOFF DO 120 I = 1,NS ICI = IC + I 120 X(ICI) = BETA2(I)*X(ICI+NS) DO 130 JX = 2,MX-1 IC = IYOFF + NS*(JX-1) DO 125 I = 1,NS ICI = IC + I 125 X(ICI) = BETA(I)*X(ICI+NS) 130 CONTINUE JX = MX IC = IYOFF + NS*(JX-1) DO 135 I = 1,NS ICI = IC + I 135 X(ICI) = 0.0E0 C----------------------------------------------------------------------- C Calculate (I - (D-inverse)*L)-inverse * X. C----------------------------------------------------------------------- 160 CONTINUE JY = 1 DO 175 JX = 2,MX-1 IC = NS*(JX-1) DO 170 I = 1,NS ICI = IC + I 170 X(ICI) = X(ICI) + BETA(I)*X(ICI-NS) 175 CONTINUE JX = MX IC = NS*(JX-1) DO 180 I = 1,NS ICI = IC + I 180 X(ICI) = X(ICI) + BETA2(I)*X(ICI-NS) DO 210 JY = 2,MY-1 IYOFF = MXNS*(JY-1) JX = 1 IC = IYOFF DO 185 I = 1,NS ICI = IC + I 185 X(ICI) = X(ICI) + GAMMA(I)*X(ICI-MXNS) DO 200 JX = 2,MX-1 IC = IYOFF + NS*(JX-1) DO 195 I = 1,NS ICI = IC + I X(ICI) = (X(ICI) + BETA(I)*X(ICI-NS)) 1 + GAMMA(I)*X(ICI-MXNS) 195 CONTINUE 200 CONTINUE JX = MX IC = IYOFF + NS*(JX-1) DO 205 I = 1,NS ICI = IC + I X(ICI) = (X(ICI) + BETA2(I)*X(ICI-NS)) 1 + GAMMA(I)*X(ICI-MXNS) 205 CONTINUE 210 CONTINUE JY = MY IYOFF = MXNS*(JY-1) JX = 1 IC = IYOFF DO 215 I = 1,NS ICI = IC + I 215 X(ICI) = X(ICI) + GAMMA2(I)*X(ICI-MXNS) DO 225 JX = 2,MX-1 IC = IYOFF + NS*(JX-1) DO 220 I = 1,NS ICI = IC + I X(ICI) = (X(ICI) + BETA(I)*X(ICI-NS)) 1 + GAMMA2(I)*X(ICI-MXNS) 220 CONTINUE 225 CONTINUE JX = MX IC = IYOFF + NS*(JX-1) DO 230 I = 1,NS ICI = IC + I X(ICI) = (X(ICI) + BETA2(I)*X(ICI-NS)) 1 + GAMMA2(I)*X(ICI-MXNS) 230 CONTINUE C----------------------------------------------------------------------- C Add increment X to Z. C----------------------------------------------------------------------- DO 300 I = 1,N 300 Z(I) = Z(I) + X(I) C IF (ITER .LT. ITMAX) GO TO 70 RETURN C------------ End of Subroutine GS ----------------------------------- END Demonstration program for SVODPK package Food web problem with NS species, NS = 6 Predator-prey interaction and diffusion on a 2-D square Matrix parameters: a = 0.1000E+01 e = 0.1000E+05 g = 0.5000E-06 b parameter = 0.1000E+01 Diffusion coefficients: Dprey = 0.1000E+01 Dpred = 0.5000E+00 Rate parameter alpha = 0.1000E+01 Mesh dimensions (mx,my) = 6 6 Total system size is NEQ = 216 Preconditioning uses interaction-only block-diagonal matrix with block-grouping, and G-S iterations Number of diagonal block groups = ngrp = 4 (ngx by ngy, ngx = 2 ngy = 2 ) G-S preconditioner uses itmax iterations, itmax = 5 Tolerance parameters: RTOL = 0.10E-03 ATOL = 0.10E-03 ************************************************************************************* Solution with MF = 10 t NSTEP NFE NNI NLI NPE NQ H AVDIM NCF rate LCF rate 0.10E-07 2 6 0 0 0 1 0.23E-07 0.00E+00 0.00E+00 0.00E+00 0.10E-06 4 11 0 0 0 2 0.67E-07 0.00E+00 0.00E+00 0.00E+00 0.10E-05 8 20 0 0 0 3 0.50E-06 0.00E+00 0.00E+00 0.00E+00 0.10E-04 18 40 0 0 0 4 0.91E-06 0.00E+00 0.00E+00 0.00E+00 0.10E-03 85 153 0 0 0 3 0.14E-05 0.00E+00 0.15E-01 0.00E+00 0.10E-02 491 921 0 0 0 1 0.33E-05 0.00E+00 0.35E+00 0.00E+00 Final statistics for this run: RWORK size = 3476 IWORK size = 30 Number of time steps = 491 Number of f evaluations = 921 Number of preconditioner evals. = 0 Number of preconditioner solves = 0 Number of nonlinear iterations = 0 Number of linear iterations = 0 Average subspace dimension = 0.0000 145 nonlinear conv. failures, 0 linear conv. failures, 2 error test failures ************************************************************************************* Solution with MF = 21 t NSTEP NFE NNI NLI NPE NQ H AVDIM NCF rate LCF rate 0.10E-07 2 7 3 1 2 1 0.23E-07 0.33E+00 0.00E+00 0.00E+00 0.10E-06 4 15 8 4 4 2 0.60E-07 0.60E+00 0.00E+00 0.00E+00 0.10E-05 11 29 15 11 7 3 0.27E-06 0.10E+01 0.00E+00 0.00E+00 0.10E-04 28 65 34 28 8 4 0.57E-06 0.89E+00 0.00E+00 0.00E+00 0.10E-03 91 202 103 96 15 4 0.72E-05 0.99E+00 0.00E+00 0.00E+00 0.10E-02 102 234 115 116 21 1 0.35E-03 0.17E+01 0.00E+00 0.00E+00 0.10E-01 109 258 124 131 23 2 0.20E-02 0.17E+01 0.00E+00 0.00E+00 0.10E+00 126 322 143 176 26 4 0.13E-01 0.24E+01 0.00E+00 0.00E+00 0.10E+01 166 544 184 357 29 4 0.34E-01 0.44E+01 0.00E+00 0.00E+00 0.20E+01 184 662 204 455 31 3 0.85E-01 0.49E+01 0.00E+00 0.90E+00 SVODPK- Warning. Poor iterative algorithm performance at T = R1. Average no. of linear iterations = R2 In above, R1 = 0.2901315212250E+01 R2 = 0.5000000000000E+01 SVODPK- Warning. Poor iterative algorithm performance at T = R1. Linear convergence failure rate = R2 In above, R1 = 0.2901315212250E+01 R2 = 0.1000000000000E+01 SVODPK- Warning. Poor iterative algorithm performance at T = R1. Average no. of linear iterations = R2 In above, R1 = 0.2985995292664E+01 R2 = 0.5000000000000E+01 SVODPK- Warning. Poor iterative algorithm performance at T = R1. Linear convergence failure rate = R2 In above, R1 = 0.2985995292664E+01 R2 = 0.1000000000000E+01 0.30E+01 196 782 224 555 31 3 0.85E-01 0.50E+01 0.00E+00 0.10E+01 0.40E+01 205 866 238 625 33 3 0.14E+00 0.50E+01 0.00E+00 0.10E+01 0.50E+01 212 914 246 665 33 3 0.14E+00 0.50E+01 0.00E+00 0.10E+01 0.60E+01 221 998 260 735 40 2 0.30E+00 0.50E+01 0.22E+00 0.79E+00 0.70E+01 224 1016 263 750 40 2 0.30E+00 0.50E+01 0.00E+00 0.10E+01 0.80E+01 227 1034 266 765 40 2 0.30E+00 0.50E+01 0.00E+00 0.10E+01 0.90E+01 230 1058 270 785 41 1 0.46E+00 0.50E+01 0.00E+00 0.10E+01 0.10E+02 232 1100 277 820 42 1 0.46E+00 0.50E+01 0.00E+00 0.10E+01 Final statistics for this run: RWORK size = 3877 IWORK size = 54 Number of time steps = 232 Number of f evaluations = 1100 Number of preconditioner evals. = 42 Number of preconditioner solves = 2168 Number of nonlinear iterations = 277 Number of linear iterations = 820 Average subspace dimension = 2.9603 2 nonlinear conv. failures, 88 linear conv. failures, 1 error test failures ************************************************************************************* Solution with MF = 29 t NSTEP NFE NNI NLI NPE NQ H AVDIM NCF rate LCF rate 0.10E-07 2 6 3 0 2 1 0.23E-07 0.00E+00 0.00E+00 0.00E+00 0.10E-06 4 11 8 0 4 2 0.60E-07 0.00E+00 0.00E+00 0.00E+00 0.10E-05 11 18 15 0 7 3 0.27E-06 0.00E+00 0.00E+00 0.00E+00 0.10E-04 28 37 34 0 8 4 0.56E-06 0.00E+00 0.00E+00 0.00E+00 0.10E-03 91 105 102 0 15 4 0.64E-05 0.00E+00 0.00E+00 0.00E+00 0.10E-02 104 120 117 0 21 1 0.36E-03 0.00E+00 0.00E+00 0.00E+00 0.10E-01 111 129 126 0 23 3 0.28E-02 0.00E+00 0.00E+00 0.00E+00 0.10E+00 126 147 144 0 26 4 0.12E-01 0.00E+00 0.00E+00 0.00E+00 0.10E+01 246 324 321 0 42 3 0.83E-02 0.00E+00 0.00E+00 0.00E+00 0.20E+01 367 504 501 0 48 3 0.83E-02 0.00E+00 0.00E+00 0.00E+00 0.30E+01 487 673 670 0 54 3 0.83E-02 0.00E+00 0.00E+00 0.00E+00 0.40E+01 608 845 842 0 60 3 0.83E-02 0.00E+00 0.00E+00 0.00E+00 0.50E+01 729 1017 1014 0 66 3 0.83E-02 0.00E+00 0.00E+00 0.00E+00 0.60E+01 849 1188 1185 0 72 3 0.83E-02 0.00E+00 0.00E+00 0.00E+00 0.70E+01 970 1362 1359 0 78 3 0.83E-02 0.00E+00 0.00E+00 0.00E+00 0.80E+01 1090 1533 1530 0 84 3 0.83E-02 0.00E+00 0.00E+00 0.00E+00 0.90E+01 1211 1707 1704 0 90 3 0.83E-02 0.00E+00 0.00E+00 0.00E+00 0.10E+02 1332 1879 1876 0 96 3 0.83E-02 0.00E+00 0.00E+00 0.00E+00 Final statistics for this run: RWORK size = 2756 IWORK size = 54 Number of time steps = 1332 Number of f evaluations = 1879 Number of preconditioner evals. = 96 Number of preconditioner solves = 3734 Number of nonlinear iterations = 1876 Number of linear iterations = 0 Average subspace dimension = 0.0000 0 nonlinear conv. failures, 0 linear conv. failures, 8 error test failures -------------------------------------------------------------------------------- At time t = 0.00000000E+00 -------------------------------------------------------------------------------- The species c( 1) values are 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.1678 10.3775 10.3775 10.1678 10.0000 10.0000 10.3775 10.8493 10.8493 10.3775 10.0000 10.0000 10.3775 10.8493 10.8493 10.3775 10.0000 10.0000 10.1678 10.3775 10.3775 10.1678 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 -------------------------------------------------------------------------------- The species c( 2) values are 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.3355 10.7550 10.7550 10.3355 10.0000 10.0000 10.7550 11.6987 11.6987 10.7550 10.0000 10.0000 10.7550 11.6987 11.6987 10.7550 10.0000 10.0000 10.3355 10.7550 10.7550 10.3355 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 -------------------------------------------------------------------------------- The species c( 3) values are 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.5033 11.1325 11.1325 10.5033 10.0000 10.0000 11.1325 12.5480 12.5480 11.1325 10.0000 10.0000 11.1325 12.5480 12.5480 11.1325 10.0000 10.0000 10.5033 11.1325 11.1325 10.5033 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 -------------------------------------------------------------------------------- The species c( 4) values are 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.6711 11.5099 11.5099 10.6711 10.0000 10.0000 11.5099 13.3974 13.3974 11.5099 10.0000 10.0000 11.5099 13.3974 13.3974 11.5099 10.0000 10.0000 10.6711 11.5099 11.5099 10.6711 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 -------------------------------------------------------------------------------- The species c( 5) values are 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.8389 11.8874 11.8874 10.8389 10.0000 10.0000 11.8874 14.2467 14.2467 11.8874 10.0000 10.0000 11.8874 14.2467 14.2467 11.8874 10.0000 10.0000 10.8389 11.8874 11.8874 10.8389 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 -------------------------------------------------------------------------------- The species c( 6) values are 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 11.0066 12.2649 12.2649 11.0066 10.0000 10.0000 12.2649 15.0961 15.0961 12.2649 10.0000 10.0000 12.2649 15.0961 15.0961 12.2649 10.0000 10.0000 11.0066 12.2649 12.2649 11.0066 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.10000000E-05 -------------------------------------------------------------------------------- The species c( 1) values are 9.99991 9.99992 9.99993 9.99994 9.99992 9.99992 9.99992 10.1677 10.3774 10.3774 10.1677 9.99992 9.99993 10.3774 10.8492 10.8492 10.3774 9.99994 9.99993 10.3774 10.8492 10.8492 10.3774 9.99993 9.99992 10.1677 10.3774 10.3774 10.1677 9.99992 9.99991 9.99992 9.99993 9.99993 9.99992 9.99991 -------------------------------------------------------------------------------- The species c( 2) values are 9.99991 9.99993 9.99995 9.99995 9.99994 9.99992 9.99993 10.3355 10.7549 10.7549 10.3355 9.99994 9.99995 10.7549 11.6985 11.6985 10.7549 9.99995 9.99995 10.7549 11.6985 11.6985 10.7549 9.99995 9.99993 10.3355 10.7549 10.7549 10.3355 9.99993 9.99991 9.99993 9.99995 9.99995 9.99993 9.99991 -------------------------------------------------------------------------------- The species c( 3) values are 9.99991 9.99994 9.99997 9.99997 9.99994 9.99992 9.99994 10.5032 11.1323 11.1323 10.5032 9.99994 9.99997 11.1323 12.5478 12.5478 11.1323 9.99997 9.99997 11.1323 12.5478 12.5478 11.1323 9.99997 9.99994 10.5032 11.1323 11.1323 10.5032 9.99994 9.99991 9.99994 9.99997 9.99997 9.99994 9.99991 -------------------------------------------------------------------------------- The species c( 4) values are 13.5005 13.5005 13.5006 13.5006 13.5005 13.5005 13.5005 14.5524 15.8954 15.8954 14.5524 13.5005 13.5006 15.8954 19.0338 19.0338 15.8954 13.5006 13.5006 15.8954 19.0338 19.0338 15.8954 13.5006 13.5005 14.5524 15.8954 15.8954 14.5524 13.5005 13.5005 13.5005 13.5006 13.5006 13.5005 13.5005 -------------------------------------------------------------------------------- The species c( 5) values are 13.5005 13.5005 13.5006 13.5006 13.5005 13.5005 13.5005 14.7812 16.4167 16.4167 14.7812 13.5005 13.5006 16.4167 20.2405 20.2405 16.4167 13.5006 13.5006 16.4167 20.2405 20.2405 16.4167 13.5006 13.5005 14.7812 16.4167 16.4167 14.7812 13.5005 13.5005 13.5005 13.5006 13.5006 13.5005 13.5005 -------------------------------------------------------------------------------- The species c( 6) values are 13.5005 13.5005 13.5006 13.5006 13.5005 13.5005 13.5005 15.0100 16.9380 16.9380 15.0100 13.5005 13.5006 16.9380 21.4471 21.4471 16.9380 13.5006 13.5006 16.9380 21.4471 21.4471 16.9380 13.5006 13.5005 15.0100 16.9380 16.9380 15.0100 13.5005 13.5005 13.5005 13.5006 13.5006 13.5005 13.5005 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.99999993E-03 -------------------------------------------------------------------------------- The species c( 1) values are 9.90731 9.91677 9.92833 9.93029 9.92263 9.91696 9.91486 10.0748 10.2771 10.2787 10.0797 9.92263 9.92447 10.2751 10.7186 10.7198 10.2787 9.93029 9.92446 10.2747 10.7177 10.7186 10.2771 9.92833 9.91483 10.0736 10.2747 10.2751 10.0748 9.91677 9.90725 9.91483 9.92446 9.92447 9.91486 9.90731 -------------------------------------------------------------------------------- The species c( 2) values are 9.90777 9.92482 9.94599 9.94798 9.93068 9.91743 9.92291 10.2414 10.6443 10.6460 10.2463 9.93068 9.94212 10.6422 11.5274 11.5288 10.6460 9.94798 9.94211 10.6418 11.5265 11.5274 10.6443 9.94599 9.92287 10.2402 10.6418 10.6422 10.2414 9.92482 9.90772 9.92287 9.94211 9.94212 9.92291 9.90777 -------------------------------------------------------------------------------- The species c( 3) values are 9.90823 9.93286 9.96367 9.96563 9.93872 9.91789 9.93095 10.4080 11.0113 11.0130 10.4129 9.93872 9.95977 11.0091 12.3349 12.3364 11.0130 9.96563 9.95975 11.0087 12.3340 12.3349 11.0113 9.96367 9.93091 10.4067 11.0087 11.0091 10.4080 9.93286 9.90818 9.93091 9.95975 9.95977 9.93095 9.90823 -------------------------------------------------------------------------------- The species c( 4) values are 297241. 297752. 298386. 298444. 297926. 297529. 297695. 307251. 319336. 319386. 307396. 297926. 298270. 319273. 345820. 345861. 319386. 298444. 298270. 319260. 345794. 345820. 319336. 298386. 297694. 307214. 319260. 319273. 307251. 297752. 297240. 297694. 298270. 298270. 297695. 297241. -------------------------------------------------------------------------------- The species c( 5) values are 297241. 297752. 298386. 298444. 297926. 297529. 297695. 307250. 319336. 319386. 307396. 297926. 298270. 319273. 345820. 345861. 319386. 298444. 298270. 319260. 345794. 345820. 319336. 298386. 297694. 307214. 319260. 319273. 307251. 297752. 297240. 297694. 298270. 298270. 297695. 297241. -------------------------------------------------------------------------------- The species c( 6) values are 297241. 297752. 298386. 298444. 297926. 297529. 297695. 307251. 319336. 319386. 307396. 297926. 298270. 319273. 345820. 345861. 319386. 298444. 298270. 319260. 345794. 345820. 319336. 298386. 297694. 307214. 319260. 319273. 307250. 297752. 297240. 297694. 298270. 298270. 297695. 297241. -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.99999994E+00 -------------------------------------------------------------------------------- The species c( 1) values are 1.58849 1.59920 1.62154 1.64765 1.67033 1.68145 1.58529 1.59502 1.61546 1.63952 1.66030 1.67032 1.57754 1.58550 1.60235 1.62233 1.63949 1.64761 1.56817 1.57409 1.58704 1.60236 1.61544 1.62148 1.56044 1.56460 1.57408 1.58544 1.59499 1.59920 1.55729 1.56044 1.56817 1.57753 1.58529 1.58847 -------------------------------------------------------------------------------- The species c( 2) values are 1.59062 1.60136 1.62361 1.64978 1.67255 1.68370 1.58743 1.59713 1.61760 1.64165 1.66251 1.67256 1.57966 1.58753 1.60453 1.62448 1.64168 1.64982 1.57029 1.57621 1.58916 1.60451 1.61762 1.62366 1.56256 1.56671 1.57621 1.58758 1.59715 1.60136 1.55941 1.56256 1.57029 1.57967 1.58744 1.59062 -------------------------------------------------------------------------------- The species c( 3) values are 1.59266 1.60341 1.62574 1.65192 1.67469 1.68584 1.58947 1.59920 1.61968 1.64378 1.66463 1.67469 1.58169 1.58961 1.60657 1.62657 1.64378 1.65192 1.57231 1.57824 1.59120 1.60657 1.61968 1.62573 1.56457 1.56873 1.57824 1.58961 1.59919 1.60341 1.56141 1.56457 1.57231 1.58169 1.58947 1.59266 -------------------------------------------------------------------------------- The species c( 4) values are 47717.2 48039.0 48707.8 49492.2 50174.1 50508.0 47621.5 47912.9 48526.5 49248.3 49873.1 50174.1 47388.4 47625.8 48133.7 48732.8 49248.3 49492.2 47107.3 47284.8 47673.3 48133.7 48526.5 48707.7 46875.3 46999.9 47284.8 47625.7 47912.7 48039.0 46780.6 46875.3 47107.2 47388.4 47621.5 47717.1 -------------------------------------------------------------------------------- The species c( 5) values are 47717.2 48039.0 48707.8 49492.2 50174.1 50508.0 47621.5 47912.8 48526.5 49248.3 49873.1 50174.1 47388.4 47625.8 48133.8 48732.8 49248.3 49492.2 47107.3 47284.8 47673.3 48133.7 48526.5 48707.7 46875.3 46999.9 47284.8 47625.7 47912.7 48039.0 46780.6 46875.3 47107.2 47388.4 47621.5 47717.1 -------------------------------------------------------------------------------- The species c( 6) values are 47717.2 48039.0 48707.8 49492.2 50174.1 50508.0 47621.5 47912.8 48526.5 49248.3 49873.1 50174.1 47388.4 47625.8 48133.7 48732.8 49248.3 49492.2 47107.3 47284.8 47673.3 48133.7 48526.5 48707.7 46875.3 46999.9 47284.8 47625.7 47912.7 48039.0 46780.6 46875.3 47107.2 47388.4 47621.5 47717.1 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.40000000E+01 -------------------------------------------------------------------------------- The species c( 1) values are 1.19550 1.20387 1.22122 1.24170 1.25941 1.26809 1.19305 1.20053 1.21653 1.23531 1.25163 1.25942 1.18670 1.19281 1.20618 1.22180 1.23537 1.24164 1.17913 1.18376 1.19421 1.20612 1.21645 1.22119 1.17303 1.17635 1.18395 1.19284 1.20045 1.20377 1.17035 1.17293 1.17912 1.18667 1.19290 1.19544 -------------------------------------------------------------------------------- The species c( 2) values are 1.19529 1.20357 1.22106 1.24153 1.25937 1.26799 1.19264 1.20025 1.21627 1.23523 1.25152 1.25936 1.18652 1.19275 1.20594 1.22175 1.23516 1.24158 1.17905 1.18368 1.19365 1.20600 1.21634 1.22108 1.17275 1.17598 1.18348 1.19270 1.20032 1.20365 1.17038 1.17285 1.17905 1.18654 1.19277 1.19532 -------------------------------------------------------------------------------- The species c( 3) values are 1.19543 1.20376 1.22118 1.24165 1.25943 1.26809 1.19288 1.20043 1.21644 1.23531 1.25162 1.25943 1.18665 1.19281 1.20610 1.22182 1.23531 1.24166 1.17912 1.18376 1.19397 1.20611 1.21645 1.22118 1.17292 1.17620 1.18376 1.19282 1.20044 1.20377 1.17039 1.17292 1.17912 1.18666 1.19290 1.19544 -------------------------------------------------------------------------------- The species c( 4) values are 35861.3 36110.8 36633.0 37247.1 37780.1 38039.3 35784.8 36011.0 36491.1 37056.9 37545.9 37780.1 35597.8 35782.6 36181.1 36652.2 37056.9 37247.0 35372.0 35510.9 35817.3 36181.1 36491.0 36633.0 35186.1 35284.3 35510.9 35782.6 36011.0 36110.8 35110.3 35186.1 35372.0 35597.8 35784.8 35861.2 -------------------------------------------------------------------------------- The species c( 5) values are 35861.3 36110.8 36633.0 37247.1 37780.1 38039.3 35784.8 36011.0 36491.0 37056.9 37545.9 37780.1 35597.8 35782.6 36181.0 36652.2 37056.9 37247.0 35372.0 35510.9 35817.3 36181.1 36491.0 36633.0 35186.1 35284.3 35510.9 35782.6 36011.0 36110.8 35110.3 35186.0 35372.0 35597.8 35784.8 35861.2 -------------------------------------------------------------------------------- The species c( 6) values are 35861.3 36110.8 36633.0 37247.1 37780.1 38039.3 35784.8 36011.0 36491.0 37056.9 37545.9 37780.1 35597.8 35782.6 36181.0 36652.2 37056.9 37247.0 35372.0 35510.9 35817.3 36181.1 36491.0 36633.0 35186.1 35284.3 35510.9 35782.6 36011.0 36110.8 35110.3 35186.0 35372.0 35597.8 35784.8 35861.2 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.70000000E+01 -------------------------------------------------------------------------------- The species c( 1) values are 1.18852 1.19681 1.21414 1.23451 1.25222 1.26082 1.18599 1.19350 1.20943 1.22821 1.24444 1.25221 1.17979 1.18592 1.19915 1.21478 1.22822 1.23454 1.17230 1.17691 1.18708 1.19913 1.20948 1.21408 1.16614 1.16939 1.17691 1.18594 1.19350 1.19682 1.16364 1.16613 1.17229 1.17975 1.18608 1.18861 -------------------------------------------------------------------------------- The species c( 2) values are 1.18864 1.19685 1.21394 1.23453 1.25220 1.26080 1.18619 1.19350 1.20943 1.22821 1.24443 1.25220 1.17981 1.18588 1.19921 1.21478 1.22817 1.23449 1.17235 1.17697 1.18713 1.19916 1.20937 1.21420 1.16608 1.16932 1.17691 1.18590 1.19349 1.19679 1.16360 1.16615 1.17230 1.17981 1.18589 1.18844 -------------------------------------------------------------------------------- The species c( 3) values are 1.18841 1.19677 1.21433 1.23449 1.25218 1.26080 1.18579 1.19349 1.20942 1.22818 1.24441 1.25219 1.17976 1.18596 1.19907 1.21476 1.22821 1.23451 1.17223 1.17684 1.18700 1.19914 1.20943 1.21414 1.16615 1.16945 1.17689 1.18592 1.19350 1.19681 1.16360 1.16608 1.17229 1.17979 1.18600 1.18853 -------------------------------------------------------------------------------- The species c( 4) values are 35654.9 35903.1 36422.7 37033.6 37564.0 37821.8 35578.8 35803.8 36281.4 36844.5 37331.1 37563.9 35392.7 35576.5 35973.0 36441.8 36844.5 37033.6 35167.9 35306.2 35611.0 35973.0 36281.4 36422.7 34982.8 35080.7 35306.2 35576.5 35803.8 35903.1 34907.5 34982.8 35167.9 35392.7 35578.8 35654.9 -------------------------------------------------------------------------------- The species c( 5) values are 35654.9 35903.1 36422.7 37033.6 37564.0 37821.8 35578.8 35803.8 36281.4 36844.5 37331.1 37564.0 35392.7 35576.5 35973.0 36441.8 36844.5 37033.6 35167.9 35306.2 35611.0 35973.0 36281.4 36422.7 34982.8 35080.7 35306.2 35576.5 35803.8 35903.1 34907.4 34982.8 35167.9 35392.7 35578.8 35654.9 -------------------------------------------------------------------------------- The species c( 6) values are 35654.9 35903.1 36422.7 37033.6 37564.0 37821.8 35578.8 35803.8 36281.4 36844.5 37331.1 37564.0 35392.7 35576.5 35973.0 36441.8 36844.5 37033.6 35167.9 35306.2 35611.0 35973.0 36281.4 36422.7 34982.8 35080.7 35306.2 35576.5 35803.8 35903.1 34907.4 34982.8 35167.9 35392.7 35578.8 35654.9 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.10000000E+02 -------------------------------------------------------------------------------- The species c( 1) values are 1.18838 1.19666 1.21399 1.23436 1.25207 1.26065 1.18584 1.19335 1.20928 1.22805 1.24429 1.25205 1.17964 1.18577 1.19900 1.21462 1.22804 1.23434 1.17213 1.17675 1.18695 1.19900 1.20928 1.21399 1.16593 1.16924 1.17677 1.18578 1.19336 1.19667 1.16341 1.16597 1.17217 1.17964 1.18585 1.18839 -------------------------------------------------------------------------------- The species c( 2) values are 1.18839 1.19668 1.21400 1.23435 1.25203 1.26066 1.18586 1.19336 1.20929 1.22804 1.24427 1.25206 1.17965 1.18579 1.19901 1.21462 1.22808 1.23442 1.17218 1.17678 1.18687 1.19900 1.20928 1.21399 1.16605 1.16928 1.17675 1.18579 1.19335 1.19666 1.16364 1.16600 1.17215 1.17966 1.18585 1.18838 -------------------------------------------------------------------------------- The species c( 3) values are 1.18838 1.19666 1.21399 1.23438 1.25205 1.26066 1.18584 1.19335 1.20928 1.22807 1.24427 1.25204 1.17964 1.18577 1.19899 1.21464 1.22805 1.23434 1.17214 1.17676 1.18696 1.19900 1.20928 1.21399 1.16596 1.16923 1.17677 1.18577 1.19335 1.19667 1.16335 1.16597 1.17215 1.17963 1.18584 1.18838 -------------------------------------------------------------------------------- The species c( 4) values are 35650.6 35898.8 36418.3 37029.2 37559.5 37817.3 35574.5 35799.5 36277.1 36840.1 37326.6 37559.5 35388.4 35572.3 35968.7 36437.4 36840.1 37029.2 35163.7 35301.9 35606.7 35968.7 36277.1 36418.4 34978.5 35076.5 35301.9 35572.3 35799.5 35898.8 34903.1 34978.5 35163.7 35388.4 35574.5 35650.6 -------------------------------------------------------------------------------- The species c( 5) values are 35650.6 35898.8 36418.3 37029.2 37559.5 37817.3 35574.5 35799.5 36277.1 36840.0 37326.6 37559.5 35388.4 35572.3 35968.7 36437.4 36840.1 37029.2 35163.7 35301.9 35606.7 35968.7 36277.1 36418.3 34978.5 35076.5 35301.9 35572.3 35799.5 35898.8 34903.1 34978.5 35163.7 35388.4 35574.5 35650.6 -------------------------------------------------------------------------------- The species c( 6) values are 35650.6 35898.8 36418.3 37029.2 37559.5 37817.3 35574.5 35799.5 36277.1 36840.0 37326.6 37559.5 35388.4 35572.3 35968.7 36437.4 36840.1 37029.2 35163.7 35301.9 35606.7 35968.7 36277.1 36418.4 34978.5 35076.5 35301.9 35572.3 35799.5 35898.8 34903.1 34978.5 35163.7 35388.4 35574.5 35650.6 --------------------------------------------------------------------------------