[2696] | 1 | MODULE EXAMPLE1 |
---|
| 2 | |
---|
| 3 | ! DEMONSTRATION PROGRAM FOR THE DVODE_F90 PACKAGE. |
---|
| 4 | |
---|
| 5 | ! The following is a simple example problem, with the coding |
---|
| 6 | ! needed for its solution by DVODE_F90. The problem is from |
---|
| 7 | ! chemical kinetics, and consists of the following three rate |
---|
| 8 | ! equations: |
---|
| 9 | ! dy1/dt = -.04d0*y1 + 1.d4*y2*y3 |
---|
| 10 | ! dy2/dt = .04d0*y1 - 1.d4*y2*y3 - 3.d7*y2**2 |
---|
| 11 | ! dy3/dt = 3.d7*y2**2 |
---|
| 12 | ! on the interval from t = 0.0d0 to t = 4.d10, with initial |
---|
| 13 | ! conditions y1 = 1.0d0, y2 = y3 = 0.0d0. The problem is stiff. |
---|
| 14 | |
---|
| 15 | ! The following coding solves this problem with DVODE_F90, |
---|
| 16 | ! using a user supplied Jacobian and printing results at |
---|
| 17 | ! t = .4, 4.,...,4.d10. It uses ITOL = 2 and ATOL much smaller |
---|
| 18 | ! for y2 than y1 or y3 because y2 has much smaller values. At |
---|
| 19 | ! the end of the run, statistical quantities of interest are |
---|
| 20 | ! printed. (See optional output in the full DVODE description |
---|
| 21 | ! below.) Output is written to the file example1.dat. |
---|
| 22 | |
---|
| 23 | CONTAINS |
---|
| 24 | |
---|
| 25 | SUBROUTINE FEX(NEQ,T,Y,YDOT) |
---|
| 26 | IMPLICIT NONE |
---|
| 27 | INTEGER NEQ |
---|
| 28 | DOUBLE PRECISION T, Y, YDOT |
---|
| 29 | DIMENSION Y(NEQ), YDOT(NEQ) |
---|
| 30 | |
---|
| 31 | YDOT(1) = -.04D0*Y(1) + 1.D4*Y(2)*Y(3) |
---|
| 32 | YDOT(3) = 3.E7*Y(2)*Y(2) |
---|
| 33 | YDOT(2) = -YDOT(1) - YDOT(3) |
---|
| 34 | RETURN |
---|
| 35 | END SUBROUTINE FEX |
---|
| 36 | |
---|
| 37 | SUBROUTINE JEX(NEQ,T,Y,ML,MU,PD,NRPD) |
---|
| 38 | IMPLICIT NONE |
---|
| 39 | INTEGER NEQ, ML, MU, NRPD |
---|
| 40 | DOUBLE PRECISION PD, T, Y |
---|
| 41 | DIMENSION Y(NEQ), PD(NRPD,NEQ) |
---|
| 42 | |
---|
| 43 | PD(1,1) = -.04D0 |
---|
| 44 | PD(1,2) = 1.D4*Y(3) |
---|
| 45 | PD(1,3) = 1.D4*Y(2) |
---|
| 46 | PD(2,1) = .04D0 |
---|
| 47 | PD(2,3) = -PD(1,3) |
---|
| 48 | PD(3,2) = 6.E7*Y(2) |
---|
| 49 | PD(2,2) = -PD(1,2) - PD(3,2) |
---|
| 50 | RETURN |
---|
| 51 | END SUBROUTINE JEX |
---|
| 52 | |
---|
| 53 | END MODULE EXAMPLE1 |
---|
| 54 | |
---|
| 55 | !****************************************************************** |
---|
| 56 | |
---|
| 57 | PROGRAM RUNEXAMPLE1 |
---|
| 58 | |
---|
| 59 | USE DVODE_F90_M |
---|
| 60 | USE EXAMPLE1 |
---|
| 61 | |
---|
| 62 | IMPLICIT NONE |
---|
| 63 | DOUBLE PRECISION ATOL, RTOL, T, TOUT, Y, RSTATS |
---|
| 64 | INTEGER NEQ, ITASK, ISTATE, ISTATS, IOUT, IERROR, I |
---|
| 65 | DIMENSION Y(3), ATOL(3), RSTATS(22), ISTATS(31) |
---|
| 66 | |
---|
| 67 | TYPE (VODE_OPTS) :: OPTIONS |
---|
| 68 | |
---|
| 69 | OPEN (UNIT=6,FILE='example1.dat') |
---|
| 70 | IERROR = 0 |
---|
| 71 | NEQ = 3 |
---|
| 72 | Y(1) = 1.0D0 |
---|
| 73 | Y(2) = 0.0D0 |
---|
| 74 | Y(3) = 0.0D0 |
---|
| 75 | T = 0.0D0 |
---|
| 76 | TOUT = 0.4D0 |
---|
| 77 | RTOL = 1.D-4 |
---|
| 78 | ATOL(1) = 1.D-8 |
---|
| 79 | ATOL(2) = 1.D-14 |
---|
| 80 | ATOL(3) = 1.D-6 |
---|
| 81 | ITASK = 1 |
---|
| 82 | ISTATE = 1 |
---|
| 83 | ! OPTIONS = SET_OPTS(DENSE_J=.TRUE.,ABSERR_VECTOR=ATOL,RELERR=RTOL, & |
---|
| 84 | ! USER_SUPPLIED_JACOBIAN=.TRUE.) |
---|
| 85 | OPTIONS = SET_NORMAL_OPTS(DENSE_J=.TRUE.,ABSERR_VECTOR=ATOL, & |
---|
| 86 | RELERR=RTOL,USER_SUPPLIED_JACOBIAN=.TRUE.) |
---|
| 87 | DO IOUT = 1, 12 |
---|
| 88 | CALL DVODE_F90(FEX,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,J_FCN=JEX) |
---|
| 89 | CALL GET_STATS(RSTATS,ISTATS) |
---|
| 90 | WRITE (6,90003) T, Y(1), Y(2), Y(3) |
---|
| 91 | DO I = 1, NEQ |
---|
| 92 | IF (Y(I)<0.0D0) IERROR = 1 |
---|
| 93 | END DO |
---|
| 94 | IF (ISTATE<0) THEN |
---|
| 95 | WRITE (6,90004) ISTATE |
---|
| 96 | STOP |
---|
| 97 | END IF |
---|
| 98 | TOUT = TOUT*10.0D0 |
---|
| 99 | END DO |
---|
| 100 | WRITE (6,90000) ISTATS(11), ISTATS(12), ISTATS(13), ISTATS(19), & |
---|
| 101 | ISTATS(20), ISTATS(21), ISTATS(22) |
---|
| 102 | IF (IERROR==1) THEN |
---|
| 103 | WRITE (6,90001) |
---|
| 104 | ELSE |
---|
| 105 | WRITE (6,90002) |
---|
| 106 | END IF |
---|
| 107 | 90000 FORMAT (/' No. steps =',I4,' No. f-s =',I4,' No. J-s =',I4, & |
---|
| 108 | ' No. LU-s =',I4/' No. nonlinear iterations =', & |
---|
| 109 | I4/' No. nonlinear convergence failures =', & |
---|
| 110 | I4/' No. error test failures =',I4/) |
---|
| 111 | 90001 FORMAT (/' An error occurred.') |
---|
| 112 | 90002 FORMAT (/' No errors occurred.') |
---|
| 113 | 90003 FORMAT (' At t =',D12.4,' y =',3D14.6) |
---|
| 114 | 90004 FORMAT (///' Error halt: ISTATE =',I3) |
---|
| 115 | STOP |
---|
| 116 | END PROGRAM RUNEXAMPLE1 |
---|