*DECK DVPKDEMO C----------------------------------------------------------------------- C Demonstration program for DVODPK. C ODE system from ns-species interaction PDE in 2 dimensions. C C This is the version of 24 April 2002. C C This version is in double 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 DVODPK 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 DGEFA and DGESL, and C the BLAS routine DAXPY. 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 DOUBLE PRECISION AA, EE, GG, BB, DPREY, DPRED, 1 AX, AY, ACOEF, BCOEF, DX, DY, ALPH, DIFF, COX, COY, 2 UROUND, SRUR DOUBLE PRECISION AVDIM, ATOL, CC, HU, RCFL, RCFN, DUMACH, 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.0D0 AY = 1.0D0 C C Call setpar to set problem parameters. CALL SETPAR C WRITE(6,10)NS 10 FORMAT(' Demonstration program for DVODPK 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 =',D12.4,' e =',D12.4, 1 ' g =',D12.4//' b parameter =',D12.4// 2 ' Diffusion coefficients: Dprey =',D12.4,' Dpred =',D12.4/) WRITE(6,16) ALPH 16 FORMAT(' Rate parameter alpha =',D12.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 = DUMACH() 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.0D-5 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 =',D10.2,' ATOL =',D10.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.0D0 TOUT = 1.0D-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 DVODPK. C DO 70 IOUT = 1,NOUT CALL DVODPK (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.0D0 RCFN = 0.0D0 RCFL = 0.0D0 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(D10.2,I5,I6,3I5,I4,2D11.2,D10.2,D12.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 =',D12.4//) GO TO 75 65 CONTINUE IF (TOUT .GT. 0.9D0) TOUT = TOUT + 1.0D0 IF (TOUT .LT. 0.9D0) TOUT = TOUT*10.0D0 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 DVODPK 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 DOUBLE PRECISION 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.0D0 EE = 1.0D4 GG = 0.5D-6 BB = 1.0D0 DPREY = 1.0D0 DPRED = 0.5D0 ALPH = 1.0D0 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.5D0 + (IG - 0.5D0)*MPER JR(NG) = 0.5D0*(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----------------------------------------------------------------------- DOUBLE PRECISION CC DIMENSION CC(*) INTEGER NS, MX, MY, MXNS INTEGER I, ICI, IOFF, IYOFF, JX, JY DOUBLE PRECISION AX,AY,ACOEF,BCOEF,DX,DY,ALPH,DIFF,COX,COY DOUBLE PRECISION 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.0D0*Y*Y*(AY-Y)*(AY-Y) IYOFF = MXNS*(JY-1) DO 10 JX = 1,MX X = (JX-1)*DX ARGX = 16.0D0*X*X*(AX-X)*(AX-X) IOFF = IYOFF + NS*(JX-1) DO 5 I = 1,NS ICI = IOFF + I CC(ICI) = 10.0D0 + 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 DOUBLE PRECISION T, C DIMENSION C(NS,MX,MY) INTEGER I, JX, JY C WRITE(LUN,10) T 10 FORMAT(/80('-')/30X,'At time t = ',D16.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 DOUBLE PRECISION T, CC, CDOT, RPAR DIMENSION CC(*), CDOT(*), RPAR(*) INTEGER NS, MX, MY, MXNS INTEGER I, IC, ICI, IDXL, IDXU, IDYL, IDYU, IYOFF, JX, JY DOUBLE PRECISION AX,AY,ACOEF,BCOEF,DX,DY,ALPH,DIFF,COX,COY DOUBLE PRECISION 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----------------------------------------------------------------------- DOUBLE PRECISION X, Y, T, C, RATE DIMENSION C(*), RATE(*) INTEGER NS, MX, MY, MXNS INTEGER I DOUBLE PRECISION AX,AY,ACOEF,BCOEF,DX,DY,ALPH,DIFF,COX,COY DOUBLE PRECISION 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.0D0 DO 15 I = 1,NS CALL DAXPY (NS, C(I), ACOEF(1,I), 1, RATE, 1) 15 CONTINUE FAC = 1.0D0 + 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 DGEFA 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 = DUMACH(). 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 DOUBLE PRECISION 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 DOUBLE PRECISION UROUND, SRUR DOUBLE PRECISION FAC, R, R0, DVNORM 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 = DVNORM (NEQ, F0, REWT) R0 = 1000.0D0*ABS(HB0)*UROUND*NEQ*FAC IF (R0 .EQ. 0.0D0) R0 = 1.0D0 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.0D0 270 IDIAG = IDIAG + (MP + 1) CALL DGEFA (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 DOUBLE PRECISION T, CC, CDOT DIMENSION CC(*), CDOT(*) INTEGER NS, MX, MY, MXNS INTEGER IBLOK, IC DOUBLE PRECISION AX,AY,ACOEF,BCOEF,DX,DY,ALPH,DIFF,COX,COY DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DGESL (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 DOUBLE PRECISION HB0, Z, X DIMENSION Z(*), X(*) INTEGER NS, MX, MY, MXNS, 1 MP, MQ, MPSQ, ITMAX INTEGER I, IC, ICI, ITER, IYOFF, JX, JY DOUBLE PRECISION AX,AY,ACOEF,BCOEF,DX,DY,ALPH,DIFF,COX,COY, 2 UROUND, SRUR DOUBLE PRECISION 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.D0/(1.D0 + 2.D0*HB0*(COX(I) + COY(I))) BETA(I) = HB0*COX(I)*ELAMDA BETA2(I) = 2.D0*BETA(I) GAMMA(I) = HB0*COY(I)*ELAMDA GAMMA2(I) = 2.D0*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.D0 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.0D0 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 DVODPK package Food web problem with NS species, NS = 6 Predator-prey interaction and diffusion on a 2-D square Matrix parameters: a = 0.1000D+01 e = 0.1000D+05 g = 0.5000D-06 b parameter = 0.1000D+01 Diffusion coefficients: Dprey = 0.1000D+01 Dpred = 0.5000D+00 Rate parameter alpha = 0.1000D+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.10D-04 ATOL = 0.10D-04 ************************************************************************************* Solution with MF = 10 t NSTEP NFE NNI NLI NPE NQ H AVDIM NCF rate LCF rate 0.10D-07 3 9 0 0 0 2 0.18D-07 0.00D+00 0.00D+00 0.00D+00 0.10D-06 5 12 0 0 0 2 0.10D-06 0.00D+00 0.00D+00 0.00D+00 0.10D-05 10 25 0 0 0 4 0.38D-06 0.00D+00 0.00D+00 0.00D+00 0.10D-04 25 41 0 0 0 5 0.65D-06 0.00D+00 0.00D+00 0.00D+00 0.10D-03 103 174 0 0 0 2 0.24D-05 0.00D+00 0.38D-01 0.00D+00 0.10D-02 490 919 0 0 0 2 0.13D-05 0.00D+00 0.36D+00 0.00D+00 Final statistics for this run: RWORK size = 3476 IWORK size = 30 Number of time steps = 490 Number of f evaluations = 919 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 142 nonlinear conv. failures, 0 linear conv. failures, 4 error test failures ************************************************************************************* Solution with MF = 21 t NSTEP NFE NNI NLI NPE NQ H AVDIM NCF rate LCF rate 0.10D-07 3 12 6 3 4 2 0.18D-07 0.50D+00 0.00D+00 0.00D+00 0.10D-06 7 20 10 7 6 2 0.48D-07 0.10D+01 0.00D+00 0.00D+00 0.10D-05 15 38 20 15 7 3 0.15D-06 0.80D+00 0.00D+00 0.00D+00 0.10D-04 34 80 43 34 10 5 0.61D-06 0.83D+00 0.00D+00 0.00D+00 0.10D-03 118 257 134 120 18 5 0.66D-05 0.95D+00 0.00D+00 0.00D+00 0.10D-02 138 314 159 152 26 2 0.42D-03 0.13D+01 0.00D+00 0.00D+00 0.10D-01 146 335 170 162 28 3 0.16D-02 0.91D+00 0.00D+00 0.00D+00 0.10D+00 169 388 196 189 32 5 0.93D-02 0.10D+01 0.00D+00 0.00D+00 0.10D+01 206 545 243 299 37 5 0.57D-01 0.23D+01 0.00D+00 0.00D+00 0.20D+01 222 622 259 360 38 5 0.92D-01 0.38D+01 0.00D+00 0.62D-01 0.30D+01 229 667 267 397 39 5 0.14D+00 0.46D+01 0.00D+00 0.50D+00 0.40D+01 235 703 273 427 40 5 0.23D+00 0.50D+01 0.00D+00 0.50D+00 0.50D+01 240 733 278 452 40 5 0.23D+00 0.50D+01 0.00D+00 0.10D+01 0.60D+01 244 763 283 477 40 5 0.23D+00 0.50D+01 0.00D+00 0.10D+01 0.70D+01 249 793 288 502 40 5 0.23D+00 0.50D+01 0.00D+00 0.10D+01 0.80D+01 253 823 293 527 40 5 0.23D+00 0.50D+01 0.00D+00 0.10D+01 0.90D+01 258 871 301 567 41 5 0.23D+00 0.50D+01 0.00D+00 0.10D+01 0.10D+02 262 901 306 592 41 5 0.23D+00 0.50D+01 0.00D+00 0.10D+01 Final statistics for this run: RWORK size = 3877 IWORK size = 54 Number of time steps = 262 Number of f evaluations = 901 Number of preconditioner evals. = 41 Number of preconditioner solves = 1730 Number of nonlinear iterations = 306 Number of linear iterations = 592 Average subspace dimension = 1.9346 0 nonlinear conv. failures, 41 linear conv. failures, 3 error test failures ************************************************************************************* Solution with MF = 29 t NSTEP NFE NNI NLI NPE NQ H AVDIM NCF rate LCF rate 0.10D-07 3 9 6 0 4 2 0.18D-07 0.00D+00 0.00D+00 0.00D+00 0.10D-06 7 13 10 0 6 2 0.48D-07 0.00D+00 0.00D+00 0.00D+00 0.10D-05 15 23 20 0 7 3 0.15D-06 0.00D+00 0.00D+00 0.00D+00 0.10D-04 34 46 43 0 10 5 0.61D-06 0.00D+00 0.00D+00 0.00D+00 0.10D-03 117 135 132 0 17 5 0.36D-05 0.00D+00 0.00D+00 0.00D+00 0.10D-02 144 166 163 0 25 2 0.61D-03 0.00D+00 0.00D+00 0.00D+00 0.10D-01 152 177 174 0 27 3 0.18D-02 0.00D+00 0.00D+00 0.00D+00 0.10D+00 176 204 201 0 30 5 0.87D-02 0.00D+00 0.00D+00 0.00D+00 0.10D+01 391 510 507 0 59 3 0.88D-02 0.00D+00 0.00D+00 0.00D+00 0.20D+01 505 683 680 0 65 3 0.88D-02 0.00D+00 0.00D+00 0.00D+00 0.30D+01 618 850 847 0 71 3 0.88D-02 0.00D+00 0.00D+00 0.00D+00 0.40D+01 731 1016 1013 0 76 3 0.88D-02 0.00D+00 0.00D+00 0.00D+00 0.50D+01 844 1181 1178 0 82 3 0.88D-02 0.00D+00 0.00D+00 0.00D+00 0.60D+01 958 1349 1346 0 88 3 0.88D-02 0.00D+00 0.00D+00 0.00D+00 0.70D+01 1071 1515 1512 0 93 3 0.88D-02 0.00D+00 0.00D+00 0.00D+00 0.80D+01 1184 1683 1680 0 99 3 0.88D-02 0.00D+00 0.00D+00 0.00D+00 0.90D+01 1297 1849 1846 0 105 3 0.88D-02 0.00D+00 0.00D+00 0.00D+00 0.10D+02 1411 2017 2014 0 110 3 0.88D-02 0.00D+00 0.00D+00 0.00D+00 Final statistics for this run: RWORK size = 2756 IWORK size = 54 Number of time steps = 1411 Number of f evaluations = 2017 Number of preconditioner evals. = 110 Number of preconditioner solves = 3998 Number of nonlinear iterations = 2014 Number of linear iterations = 0 Average subspace dimension = 0.0000 0 nonlinear conv. failures, 0 linear conv. failures, 14 error test failures -------------------------------------------------------------------------------- At time t = 0.00000000D+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.10000000D-05 -------------------------------------------------------------------------------- The species c( 1) values are 9.99991 9.99992 9.99993 9.99993 9.99993 9.99992 9.99992 10.1677 10.3774 10.3774 10.1677 9.99993 9.99993 10.3774 10.8492 10.8492 10.3774 9.99993 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.99993 9.99992 9.99993 10.3355 10.7549 10.7549 10.3355 9.99993 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.4990 13.4990 13.4990 13.4990 13.4990 13.4989 13.4990 14.5506 15.8932 15.8932 14.5506 13.4990 13.4990 15.8932 19.0308 19.0308 15.8932 13.4990 13.4990 15.8932 19.0308 19.0308 15.8932 13.4990 13.4990 14.5506 15.8932 15.8932 14.5506 13.4990 13.4990 13.4990 13.4990 13.4990 13.4990 13.4990 -------------------------------------------------------------------------------- The species c( 5) values are 13.4990 13.4990 13.4990 13.4990 13.4990 13.4989 13.4990 14.7794 16.4145 16.4145 14.7794 13.4990 13.4990 16.4145 20.2373 20.2373 16.4145 13.4990 13.4990 16.4145 20.2373 20.2373 16.4145 13.4990 13.4990 14.7794 16.4145 16.4145 14.7794 13.4990 13.4990 13.4990 13.4990 13.4990 13.4990 13.4990 -------------------------------------------------------------------------------- The species c( 6) values are 13.4990 13.4990 13.4990 13.4990 13.4990 13.4989 13.4990 15.0082 16.9357 16.9357 15.0082 13.4990 13.4990 16.9357 21.4437 21.4437 16.9357 13.4990 13.4990 16.9357 21.4437 21.4437 16.9357 13.4990 13.4990 15.0082 16.9357 16.9357 15.0082 13.4990 13.4990 13.4990 13.4990 13.4990 13.4990 13.4990 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.10000000D-02 -------------------------------------------------------------------------------- The species c( 1) values are 9.90702 9.91664 9.92836 9.93033 9.92253 9.91674 9.91472 10.0746 10.2769 10.2785 10.0795 9.92253 9.92446 10.2748 10.7181 10.7194 10.2785 9.93033 9.92445 10.2744 10.7173 10.7181 10.2769 9.92836 9.91469 10.0734 10.2744 10.2748 10.0746 9.91664 9.90697 9.91469 9.92445 9.92446 9.91472 9.90702 -------------------------------------------------------------------------------- The species c( 2) values are 9.90741 9.92474 9.94623 9.94820 9.93064 9.91713 9.92282 10.2412 10.6440 10.6457 10.2461 9.93064 9.94232 10.6419 11.5267 11.5281 10.6457 9.94820 9.94231 10.6415 11.5258 11.5267 10.6440 9.94623 9.92279 10.2400 10.6415 10.6419 10.2412 9.92474 9.90736 9.92279 9.94231 9.94232 9.92282 9.90741 -------------------------------------------------------------------------------- The species c( 3) values are 9.90780 9.93284 9.96409 9.96606 9.93874 9.91752 9.93092 10.4078 11.0109 11.0127 10.4127 9.93874 9.96018 11.0088 12.3339 12.3354 11.0127 9.96606 9.96017 11.0083 12.3329 12.3339 11.0109 9.96409 9.93089 10.4065 11.0083 11.0088 10.4078 9.93284 9.90776 9.93089 9.96017 9.96018 9.93092 9.90780 -------------------------------------------------------------------------------- The species c( 4) values are 297231. 297749. 298393. 298451. 297925. 297520. 297692. 307244. 319327. 319378. 307390. 297925. 298276. 319264. 345799. 345840. 319378. 298451. 298276. 319252. 345771. 345799. 319327. 298393. 297691. 307208. 319252. 319264. 307244. 297749. 297229. 297691. 298276. 298276. 297692. 297231. -------------------------------------------------------------------------------- The species c( 5) values are 297231. 297749. 298393. 298451. 297925. 297520. 297692. 307244. 319327. 319378. 307390. 297925. 298276. 319264. 345799. 345840. 319378. 298451. 298276. 319252. 345771. 345799. 319327. 298393. 297691. 307208. 319252. 319264. 307244. 297749. 297229. 297691. 298276. 298276. 297692. 297231. -------------------------------------------------------------------------------- The species c( 6) values are 297231. 297749. 298393. 298451. 297925. 297520. 297692. 307244. 319327. 319378. 307390. 297925. 298276. 319264. 345799. 345840. 319378. 298451. 298276. 319252. 345771. 345799. 319327. 298393. 297691. 307208. 319252. 319264. 307244. 297749. 297229. 297691. 298276. 298276. 297692. 297231. -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.10000000D+01 -------------------------------------------------------------------------------- The species c( 1) values are 1.58859 1.59932 1.62159 1.64773 1.67044 1.68157 1.58541 1.59511 1.61556 1.63960 1.66041 1.67044 1.57764 1.58555 1.60247 1.62243 1.63960 1.64773 1.56828 1.57420 1.58714 1.60247 1.61556 1.62159 1.56056 1.56471 1.57420 1.58555 1.59511 1.59932 1.55740 1.56056 1.56828 1.57764 1.58541 1.58859 -------------------------------------------------------------------------------- The species c( 2) values are 1.59075 1.60149 1.62379 1.64995 1.67269 1.68384 1.58756 1.59728 1.61775 1.64182 1.66265 1.67269 1.57979 1.58771 1.60465 1.62463 1.64182 1.64995 1.57042 1.57634 1.58930 1.60465 1.61775 1.62379 1.56269 1.56684 1.57634 1.58771 1.59728 1.60149 1.55953 1.56269 1.57042 1.57979 1.58756 1.59075 -------------------------------------------------------------------------------- The species c( 3) values are 1.59279 1.60354 1.62587 1.65206 1.67482 1.68598 1.58960 1.59932 1.61981 1.64391 1.66477 1.67482 1.58182 1.58974 1.60670 1.62670 1.64391 1.65206 1.57244 1.57837 1.59134 1.60670 1.61981 1.62587 1.56470 1.56886 1.57837 1.58974 1.59932 1.60354 1.56154 1.56470 1.57244 1.58182 1.58960 1.59279 -------------------------------------------------------------------------------- The species c( 4) values are 47720.8 48042.7 48711.5 49496.0 50178.0 50511.9 47625.2 47916.4 48530.2 49252.1 49876.9 50178.0 47392.1 47629.4 48137.4 48736.6 49252.1 49496.0 47110.9 47288.5 47677.0 48137.4 48530.2 48711.5 46879.0 47003.5 47288.5 47629.4 47916.4 48042.7 46784.2 46879.0 47110.9 47392.1 47625.2 47720.8 -------------------------------------------------------------------------------- The species c( 5) values are 47720.8 48042.7 48711.5 49496.0 50178.0 50511.9 47625.2 47916.4 48530.2 49252.1 49876.9 50178.0 47392.1 47629.4 48137.4 48736.6 49252.1 49496.0 47110.9 47288.5 47677.0 48137.4 48530.2 48711.5 46879.0 47003.5 47288.5 47629.4 47916.4 48042.7 46784.2 46879.0 47110.9 47392.1 47625.2 47720.8 -------------------------------------------------------------------------------- The species c( 6) values are 47720.8 48042.7 48711.5 49496.0 50178.0 50511.9 47625.2 47916.4 48530.2 49252.1 49876.9 50178.0 47392.1 47629.4 48137.4 48736.6 49252.1 49496.0 47110.9 47288.5 47677.0 48137.4 48530.2 48711.5 46879.0 47003.5 47288.5 47629.4 47916.4 48042.7 46784.2 46879.0 47110.9 47392.1 47625.2 47720.8 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.40000000D+01 -------------------------------------------------------------------------------- The species c( 1) values are 1.19536 1.20369 1.22110 1.24158 1.25936 1.26801 1.19281 1.20036 1.21637 1.23524 1.25154 1.25936 1.18658 1.19274 1.20603 1.22174 1.23524 1.24158 1.17905 1.18369 1.19390 1.20603 1.21637 1.22110 1.17285 1.17613 1.18369 1.19274 1.20036 1.20369 1.17033 1.17285 1.17905 1.18658 1.19281 1.19536 -------------------------------------------------------------------------------- The species c( 2) values are 1.19539 1.20372 1.22114 1.24161 1.25939 1.26804 1.19284 1.20039 1.21640 1.23527 1.25158 1.25939 1.18661 1.19278 1.20607 1.22178 1.23527 1.24161 1.17908 1.18372 1.19393 1.20607 1.21640 1.22114 1.17288 1.17617 1.18372 1.19278 1.20039 1.20372 1.17036 1.17288 1.17908 1.18661 1.19284 1.19539 -------------------------------------------------------------------------------- The species c( 3) values are 1.19542 1.20375 1.22117 1.24165 1.25942 1.26807 1.19287 1.20042 1.21643 1.23530 1.25161 1.25942 1.18664 1.19281 1.20610 1.22181 1.23530 1.24165 1.17912 1.18375 1.19397 1.20610 1.21643 1.22117 1.17292 1.17620 1.18375 1.19281 1.20042 1.20375 1.17039 1.17292 1.17912 1.18664 1.19287 1.19542 -------------------------------------------------------------------------------- The species c( 4) values are 35860.8 36110.4 36632.6 37246.7 37779.7 38038.9 35784.4 36010.6 36490.6 37056.5 37545.5 37779.7 35597.4 35782.2 36180.7 36651.8 37056.5 37246.7 35371.6 35510.6 35816.9 36180.7 36490.6 36632.6 35185.6 35284.0 35510.6 35782.2 36010.6 36110.4 35109.9 35185.6 35371.6 35597.4 35784.4 35860.8 -------------------------------------------------------------------------------- The species c( 5) values are 35860.8 36110.4 36632.6 37246.7 37779.7 38038.9 35784.4 36010.6 36490.6 37056.5 37545.5 37779.7 35597.4 35782.2 36180.7 36651.8 37056.5 37246.7 35371.6 35510.6 35816.9 36180.7 36490.6 36632.6 35185.6 35284.0 35510.6 35782.2 36010.6 36110.4 35109.9 35185.6 35371.6 35597.4 35784.4 35860.8 -------------------------------------------------------------------------------- The species c( 6) values are 35860.8 36110.4 36632.6 37246.7 37779.7 38038.9 35784.4 36010.6 36490.6 37056.5 37545.5 37779.7 35597.4 35782.2 36180.7 36651.8 37056.5 37246.7 35371.6 35510.6 35816.9 36180.7 36490.6 36632.6 35185.6 35284.0 35510.6 35782.2 36010.6 36110.4 35109.9 35185.6 35371.6 35597.4 35784.4 35860.8 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.70000000D+01 -------------------------------------------------------------------------------- The species c( 1) values are 1.18854 1.19683 1.21416 1.23453 1.25222 1.26082 1.18600 1.19351 1.20944 1.22822 1.24445 1.25222 1.17980 1.18594 1.19916 1.21478 1.22822 1.23453 1.17231 1.17692 1.18709 1.19916 1.20944 1.21416 1.16613 1.16940 1.17692 1.18594 1.19351 1.19683 1.16362 1.16613 1.17231 1.17980 1.18600 1.18854 -------------------------------------------------------------------------------- The species c( 2) values are 1.18854 1.19683 1.21416 1.23453 1.25222 1.26083 1.18600 1.19351 1.20944 1.22822 1.24445 1.25222 1.17980 1.18593 1.19916 1.21480 1.22822 1.23453 1.17231 1.17692 1.18708 1.19916 1.20944 1.21416 1.16614 1.16940 1.17692 1.18593 1.19351 1.19683 1.16363 1.16614 1.17231 1.17980 1.18600 1.18854 -------------------------------------------------------------------------------- The species c( 3) values are 1.18854 1.19683 1.21416 1.23453 1.25222 1.26083 1.18600 1.19351 1.20944 1.22823 1.24445 1.25222 1.17980 1.18593 1.19916 1.21481 1.22823 1.23453 1.17231 1.17692 1.18708 1.19916 1.20944 1.21416 1.16615 1.16940 1.17692 1.18593 1.19351 1.19683 1.16365 1.16615 1.17231 1.17980 1.18600 1.18854 -------------------------------------------------------------------------------- The species c( 4) values are 35655.3 35903.6 36423.2 37034.2 37564.6 37822.4 35579.1 35804.2 36281.9 36845.1 37331.7 37564.6 35393.0 35576.9 35973.5 36442.5 36845.1 37034.2 35168.4 35306.6 35611.4 35973.5 36281.9 36423.2 34983.3 35081.1 35306.6 35576.9 35804.2 35903.6 34908.0 34983.3 35168.4 35393.0 35579.1 35655.3 -------------------------------------------------------------------------------- The species c( 5) values are 35655.3 35903.6 36423.2 37034.2 37564.6 37822.4 35579.1 35804.2 36281.9 36845.1 37331.7 37564.6 35393.0 35576.9 35973.5 36442.5 36845.1 37034.2 35168.4 35306.6 35611.4 35973.5 36281.9 36423.2 34983.3 35081.1 35306.6 35576.9 35804.2 35903.6 34908.0 34983.3 35168.4 35393.0 35579.1 35655.3 -------------------------------------------------------------------------------- The species c( 6) values are 35655.3 35903.6 36423.2 37034.2 37564.6 37822.4 35579.1 35804.2 36281.9 36845.1 37331.7 37564.6 35393.0 35576.9 35973.5 36442.5 36845.1 37034.2 35168.4 35306.6 35611.4 35973.5 36281.9 36423.2 34983.3 35081.1 35306.6 35576.9 35804.2 35903.6 34908.0 34983.3 35168.4 35393.0 35579.1 35655.3 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- At time t = 0.10000000D+02 -------------------------------------------------------------------------------- The species c( 1) values are 1.18839 1.19667 1.21399 1.23437 1.25204 1.26065 1.18586 1.19336 1.20928 1.22806 1.24428 1.25204 1.17966 1.18579 1.19900 1.21462 1.22806 1.23437 1.17216 1.17677 1.18693 1.19900 1.20928 1.21399 1.16598 1.16925 1.17677 1.18579 1.19336 1.19667 1.16344 1.16598 1.17216 1.17966 1.18586 1.18839 -------------------------------------------------------------------------------- The species c( 2) values are 1.18838 1.19667 1.21399 1.23437 1.25205 1.26066 1.18585 1.19336 1.20928 1.22805 1.24428 1.25205 1.17964 1.18578 1.19900 1.21462 1.22805 1.23437 1.17215 1.17677 1.18693 1.19900 1.20928 1.21399 1.16598 1.16924 1.17677 1.18578 1.19336 1.19667 1.16347 1.16598 1.17215 1.17964 1.18585 1.18838 -------------------------------------------------------------------------------- The species c( 3) values are 1.18838 1.19666 1.21399 1.23436 1.25206 1.26067 1.18583 1.19335 1.20928 1.22804 1.24427 1.25206 1.17963 1.18577 1.19900 1.21462 1.22804 1.23436 1.17214 1.17676 1.18693 1.19900 1.20928 1.21399 1.16598 1.16924 1.17676 1.18577 1.19335 1.19666 1.16350 1.16598 1.17214 1.17963 1.18583 1.18838 -------------------------------------------------------------------------------- The species c( 4) values are 35650.7 35898.9 36418.3 37029.2 37559.4 37817.4 35574.5 35799.6 36277.1 36839.9 37326.5 37559.4 35388.4 35572.4 35968.8 36437.3 36839.9 37029.2 35163.6 35301.9 35606.8 35968.8 36277.1 36418.3 34978.4 35076.4 35301.9 35572.4 35799.6 35898.9 34903.2 34978.4 35163.6 35388.4 35574.5 35650.7 -------------------------------------------------------------------------------- The species c( 5) values are 35650.7 35898.9 36418.3 37029.2 37559.4 37817.4 35574.5 35799.6 36277.1 36839.9 37326.5 37559.4 35388.4 35572.4 35968.8 36437.3 36839.9 37029.2 35163.6 35301.9 35606.8 35968.8 36277.1 36418.3 34978.4 35076.4 35301.9 35572.4 35799.6 35898.9 34903.2 34978.4 35163.6 35388.4 35574.5 35650.7 -------------------------------------------------------------------------------- The species c( 6) values are 35650.7 35898.9 36418.3 37029.2 37559.4 37817.4 35574.5 35799.6 36277.1 36839.9 37326.5 37559.4 35388.4 35572.4 35968.8 36437.3 36839.9 37029.2 35163.6 35301.9 35606.8 35968.8 36277.1 36418.3 34978.4 35076.4 35301.9 35572.4 35799.6 35898.9 34903.2 34978.4 35163.6 35388.4 35574.5 35650.7 --------------------------------------------------------------------------------