MODULE EXAMPLE1 ! DEMONSTRATION PROGRAM FOR THE DVODE_F90 PACKAGE. ! The following is a simple example problem, with the coding ! needed for its solution by DVODE_F90. The problem is from ! chemical kinetics, and consists of the following three rate ! equations: ! dy1/dt = -.04d0*y1 + 1.d4*y2*y3 ! dy2/dt = .04d0*y1 - 1.d4*y2*y3 - 3.d7*y2**2 ! dy3/dt = 3.d7*y2**2 ! on the interval from t = 0.0d0 to t = 4.d10, with initial ! conditions y1 = 1.0d0, y2 = y3 = 0.0d0. The problem is stiff. ! The following coding solves this problem with DVODE_F90, ! using a user supplied Jacobian and printing results at ! t = .4, 4.,...,4.d10. It uses ITOL = 2 and ATOL much smaller ! for y2 than y1 or y3 because y2 has much smaller values. At ! the end of the run, statistical quantities of interest are ! printed. (See optional output in the full DVODE description ! below.) Output is written to the file example1.dat. CONTAINS SUBROUTINE FEX(NEQ,T,Y,YDOT) IMPLICIT NONE INTEGER NEQ DOUBLE PRECISION T, Y, YDOT DIMENSION Y(NEQ), YDOT(NEQ) intent(in) :: NEQ, T, Y intent(out) :: YDOT ! COMMON /CV/EXI,C0,C1,C2,C3,C4,C5,GAM,BETA,EP & ! ,RM,GU,QQ,SM,DEL0,RMT,XM,RE,TAU REAL:: EXI,C0,C1,C2,C3,C4,GAM,BETA,GU,EP,ALM,XM,PM,PP,SM,RM,DEL0,QQ,TAU EXI=0.1d0 C0=EXI-5.d0/2.d0 C1=EXI/2.D0-1.d0/4.d0 C2=-0.5d0 C3=-0.5d0 C4=EXI-3.d0/2.d0 GAM=5.d0/3.d00 BETA=EXI/2.d0+3.D0/4.D0 GU=5.7D-3 EP=0.1D0 ALM=1.d0 XM=1.0D0 PM=1.D0 PP=PM*EP+3*XM/ALM**2 SM=ALM*PP*GU**0.5 RM=PP/EP DEL0=(1-EP**2*(SM**2/2.D0+2*(2-BETA)+GU*RM))**0.5 QQ=ALM*DEL0*(PP-PM*EP)/2.D0/GU**0.5 TAU=PP/PM/EP-1.D0 !!!!!!!!!!!!!VARIABLE DI'INTEGRATION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! YDOT(1) = 1.D0 !!!!!!!!!!!!RESISTIVITE MAGNETIQUE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! YDOT(10)=-2.D0*Y(1)*EXP(-Y(1)**2) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!VITESSE TOROIDALE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! YDOT(2)=(Y(4)*Y(3)*Y(2)+ & TAU/(1.D0+TAU)*(Y(9)*Y(6)/Y(10)-C1/BETA*Y(7)*Y(8)) & -1.D0/(1.D0+TAU)*Y(3)*Y(10)) & /(2.D0*Y(3)*(Y(1)*Y(4)+Y(5))) !!!!!!!!!!!!!!!!!DENSITE VOLUMIQUE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! YDOT(3)=((EXI-1.D0)*Y(3)*Y(4)*EP**2*SM**2*Y(3)*(Y(1)*Y(4)+Y(5)) & -Y(1)*Y(3)*& (C2*SM**2*EP**2*Y(3)*Y(4)**2-Y(3)*DEL0**2*Y(2)**2 & +Y(3)/(1.D0+EP**2*Y(1)**2)**1.5D0+EP**2*C0*Y(11) & +GU*QQ**2*EP**2*Y(7)*(C1*Y(7)-Y(1)*Y(6)/Y(10)) & +GU/BETA*(Y(9)-Y(1)/BETA*Y(8))* & (-RM*EP**2/Y(10))*(BETA*Y(4)*Y(9)-Y(8)*(Y(1)*Y(4)+Y(5)))) & -Y(3)* & (C3*SM**2*EP**2*Y(3)*Y(4)*Y(5)-Y(1)*Y(3)/ & (1.D0+EP**2*Y(1)**2)**1.5D0 & -GU*QQ**2*Y(7)*Y(6)/Y(10)-GU/BETA**2/EP**2*Y(8)* & (-RM*EP**2/Y(10))*(BETA*Y(4)*Y(9)-Y(8)*(Y(1)*Y(4)+Y(5)))) & -1.D0/BETA*Y(3)*Y(9)**(-1.D0-1.D0/BETA)*Y(8)* & Y(3)*(1.D0+Y(1)**2*EP**2))& /(Y(3)*(SM**2*EP**2*(Y(1)*Y(4)+Y(5))**2-Y(9)**(-1.D0/BETA)* & (1.D0+EP**2*Y(1)**2))) !!!!!!!!!!!!!!!!!!PRESSION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! YDOT(11)=-1.D0/BETA*Y(3)*Y(9)**(-1.D0-1.D0/BETA)*Y(8)+ & Y(9)**(-1.D0/BETA)*YDOT(3) !!!!!!!!!!!!!!!!!EQUILIBRE RADIALE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! YDOT(4)=(C2*SM**2*EP**2*Y(3)*Y(4)**2 & -Y(3)*DEL0**2*Y(2)**2 & +Y(3)/(1.D0+(EP*Y(1))**2)**1.5D0 & +EP**2*C0*Y(11) & +GU*QQ**2*EP**2*Y(7)*(C1*Y(7)-Y(1)*Y(6)/Y(10)) & +GU/BETA*(Y(9)-Y(1)*Y(8)/BETA)* & (-RM*EP**2/Y(10))*(BETA*Y(4)*Y(9)-Y(8)*(Y(1)*Y(4)+Y(5))) & -EP**2*Y(1)*YDOT(11)) & /(EP**2*SM**2*Y(3)*(Y(1)*Y(4)+Y(5))) !!!!!!!!!!!!!!!!!!EQUILIBRE VERTICALE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! YDOT(5)=(C3*EP**2*SM**2*Y(3)*Y(4)*Y(5)- & Y(1)*Y(3)/(1.D0+(EP*Y(1))**2)**1.5D0 & -GU*QQ**2*Y(7)*Y(6)/Y(10)- & GU*Y(8)/BETA**2/EP**2* & (-RM*EP**2/Y(10))*(BETA*Y(4)*Y(9)-Y(8)*(Y(1)*Y(4)+Y(5))) & -YDOT(11)) & /(EP**2*SM**2*Y(3)*(Y(1)*Y(4)+Y(5))) !!!!!!!!!!!!!!!!!!INDUCTION MAGNETIQUE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! YDOT(6)= (EP**2*Y(1)*((C1-1.D0)*Y(6)+C1*Y(7)*YDOT(10)) & -EP**2*(C1*Y(7)*Y(10)-Y(1)*Y(6))*(C1-1.5D0)& -XM*RM*DEL0/QQ/SM*(1.5D0/BETA*Y(8)*Y(2)+Y(9)*YDOT(2)) & +XM*RM*EP**2*BETA*Y(7)*Y(4) & +XM*RM*EP**2*(Y(5)+Y(1)*Y(4))* & (Y(6)/Y(10)-Y(7)*YDOT(3)/Y(3))) & /(1.D0+EP**2*Y(1)**2) YDOT(7)=Y(6)/Y(10) !!!!!!!!!!!!!!!!!FLUX MAGNETIQUE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! YDOT(8)=(EP**2*(BETA*(2.D0-BETA)*Y(9)+(2*BETA-3.D0)*Y(1)*Y(8)) & -RM*EP**2/Y(10)*(BETA*Y(4)*Y(9)-Y(8)*(Y(1)*Y(4)+Y(5)))) & /(1.D0+(EP*Y(1))**2) YDOT(9)=Y(8) RETURN END SUBROUTINE FEX SUBROUTINE JEX(NEQ,T,Y,ML,MU,PD,NRPD) IMPLICIT NONE INTEGER NEQ, ML, MU, NRPD DOUBLE PRECISION PD, T, Y DIMENSION Y(NEQ), PD(NRPD,NEQ) ! PD(1,1) = -.04D0 ! PD(1,2) = 1.D4*Y(3) ! PD(1,3) = 1.D4*Y(2) ! PD(2,1) = .04D0 ! PD(2,3) = -PD(1,3) ! PD(3,2) = 6.E7*Y(2) ! PD(2,2) = -PD(1,2) - PD(3,2) RETURN END SUBROUTINE JEX END MODULE EXAMPLE1 !****************************************************************** PROGRAM RUNEXAMPLE1 USE DVODE_F90_M USE EXAMPLE1 IMPLICIT NONE DOUBLE PRECISION ATOL, RTOL, T, TOUT, Y, RSTATS INTEGER NEQ, ITASK, ISTATE, ISTATS, IOUT, IERROR, I DIMENSION Y(11), ATOL(11), RSTATS(22), ISTATS(31) TYPE (VODE_OPTS) :: OPTIONS ! COMMON /CV/EXI,C0,C1,C2,C3,C4,C5,GAM,BETA,EP & ! ,RM,GU,QQ,SM,DEL0,RMT,XM,RE,TAU OPEN (UNIT=6,FILE='example1.dat') IERROR = 0 NEQ = 11 Y(1) = 0.0001D0 Y(2) = 1.0D0 Y(3) = 1.0D0 Y(4) = 1.0D0 Y(5) = 0.0D0 Y(6) = -1.D0 Y(7) = 0.0D0 Y(8) = 0.0D0 Y(9) = 1.0D0 Y(10) = 1.0D0 Y(11) = 1.0D0 T = 0.0D0 TOUT = 0.4D0 RTOL = 1.D-4 ATOL(1) = 1.D-8 ATOL(2) = 1.D-14 ATOL(3) = 1.D-6 ATOL(4) = 1.D-6 ATOL(5) = 1.D-6 ATOL(6) = 1.D-6 ATOL(7) = 1.D-6 ATOL(8) = 1.D-6 ATOL(9) = 1.D-6 ATOL(10) = 1.D-6 ATOL(11) = 1.D-6 ITASK = 1 ISTATE = 1 OPTIONS = SET_NORMAL_OPTS(DENSE_J=.TRUE.,ABSERR_VECTOR=ATOL, & RELERR=RTOL,USER_SUPPLIED_JACOBIAN=.FALSE.) DO IOUT = 1, 12 CALL DVODE_F90(FEX,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,J_FCN=JEX) CALL GET_STATS(RSTATS,ISTATS) WRITE (6,90003) T, Y(1), Y(2), Y(3) DO I = 1, NEQ IF (Y(I)<0.0D0) IERROR = 1 END DO IF (ISTATE<0) THEN WRITE (6,90004) ISTATE STOP END IF TOUT = TOUT*10.0D0 END DO WRITE (6,90000) ISTATS(11), ISTATS(12), ISTATS(13), ISTATS(19), & ISTATS(20), ISTATS(21), ISTATS(22) IF (IERROR==1) THEN WRITE (6,90001) ELSE WRITE (6,90002) END IF 90000 FORMAT (/' No. steps =',I4,' No. f-s =',I4,' No. J-s =',I4, & ' No. LU-s =',I4/' No. nonlinear iterations =', & I4/' No. nonlinear convergence failures =', & I4/' No. error test failures =',I4/) 90001 FORMAT (/' An error occurred.') 90002 FORMAT (/' No errors occurred.') 90003 FORMAT (' At t =',D12.4,' y =',3D14.6) 90004 FORMAT (///' Error halt: ISTATE =',I3) STOP END PROGRAM RUNEXAMPLE1