Fortran90 Diffusion Program
PROGRAM DIFF2
IMPLICIT NONE
!
! THIS PROGRAM USES THE LAPLACE EQUATION TO SOLVE FOR
! DIFFUSION. THE PROGRAM USES A TECHNIQUE CALLED SUCCESSIVE
! OVER RELAXATION (SOR). THE PROGRAM CAN USE ANY RECTANGULAR
! GRID WITH A BOUNDARY ON THAT GRID. THE VALUES AT THE
! BOUNDARY MUST BE GIVEN ALONG WITH GUESSES FOR INTERIOR
! POINT VALUES. OTHER INPUT PARAMETERS ARE EXPLAINED
! THROUGHOUT THE PROGRAM.
!
! UU = ARRAY IN (X,Y) FORM
REAL, DIMENSION(51,51) :: UU
!
! ILHS = INTEGER LEFT HAND SIDE - LEFT HAND SIDE TO INTEGRATE
! IRHS = INTEGER RIGHT HAND SIDE - RIGHT HAND SIDE TO INTEGRATE
! A, B, C, D, E, G, H = WEIGHTED AVERAGES
! U = ARRAY OF ALL POINTS
REAL, DIMENSION(2601) :: A, B, C, D, E, F,G, H, U
INTEGER, DIMENSION(49) :: ILHS, IRHS
!
INTEGER :: N, LVERT, LHORZ, TT, ITMAX, ITER, NM, IC, M, ICP1&
&, ITIME, K,IBM, IB, ICOUNT, II, JJ, IZ, IA, JB
REAL :: W, SIGMA, THETA, PAR, ERRSOR, DH, DT, TEND, TD, DNRM
!
! OPEN OUTPUT FILES
!
OPEN(6,FILE='DIFF2.OUT')
OPEN(7,FILE='TIMES2.DAT')
OPEN(8,FILE='CONCEN2.ASCII')
!
! N = THE NUMBER OF INTERVALS IN A UNIT LENGTH
N=25
! LVERT= THE NUMBER OF UNIT LENGTHS ON THE VERTICAL
! SIDE OF THE GRID - IN THIS CASE 1 MM LONG
LVERT=2
! LHORIZ = THE NUMBER OF UNIT LENGTHS ON THE HORIZONTAL
! SIDE OF THE GRID
LHORZ=2
! W = CONVERGENCE FACTOR FOR THE SUCCESSIVE
! OVER-RELAXATION METHOD
W=1.85
! SIGMA = THE CONSTANT IN THE LAPLACE EQUATION -
! IN THIS CASE THE DIFFUSION CONSTANT
SIGMA=1.6/(10**5)
! THETA = A PARAMETER SUCH THAT 1/2 < THETA < 1 -
! SO THAT THE VALUES WILL ALWAYS CONVERGE
THETA=0.75
! PAR = (A CHANGE IN TIME)/((1/N)**2)
PAR=15875.0
! TT = THE TOTAL NUMBER OF TIME STEPS
TT=1000
! ITMAX = THE LARGEST NUMBER OF ITERATIONS FOR CONVERGENCE.
ITMAX=500
! ERRSOR = THE SOR ERROR TO EXPECT IN THE FINAL SOLUTIONS.
ERRSOR=0.001
!
! DELTA LENGTH
DH=1./N
! DELTA TIME
DT=PAR*(DH**2)
! TOTAL TIME
TEND=TT*DT
NM=LVERT*N-1
IC=LHORZ*N+1
M=IC*(NM+2)
ICP1=IC+1
WRITE(6,19)DH,DT,TEND,NM,IC,M
19 FORMAT(3E10.3,3I10)
!
! ASSIGNING INTERIOR POINT VALUES
DO JJ=1,M
U(JJ)=0.25
END DO
!
! LOCATING BOUNDARY AND ASSIGNING VALUES.
!
DO K=1, NM
ILHS(K)=(K-1)*IC+IC+2
IRHS(K)=(K-1)*IC+2*IC-1
END DO
!
TD=0.0
IBM=(IC-1)/2+1
DO IB=2,50
U(IB)=1.20
END DO
!
! ESTABLISH COEFF.S FOR INTERATION
!
CALL ESTAB(N, LHORZ, LVERT ,SIGMA,THETA, DT, DH&
&, A, B, C, D, E, F,G, H)
!
! BEGIN TIME LOOP
!
ITIME=1
!
! WRITE CONCENTRATIONS FOR EACH TIME STEP
! PUT U INTO ARRAY OF (X,Y)
28 ICOUNT=1
II=1
JJ=1
DO IZ=1,M
IF(ICOUNT.EQ.ICP1)II=1
IF(ICOUNT.EQ.ICP1)JJ=JJ+1
IF(ICOUNT.EQ.ICP1)ICOUNT=1
UU(II,JJ)=U(IZ)
ICOUNT=ICOUNT+1
II=II+1
END DO
WRITE(8,95)((UU(IA,JB),IA=1,IC),JB=1,IC)
WRITE(6,203)ITIME,TD,UU(IBM,12)
WRITE(7,204)ITIME,TD,UU(IBM,12)
!
ITIME=ITIME+1
TD=TD+DT
CALL SLEEP(N,LHORZ,LVERT,W,U,ERRSOR,ITMAX,DNRM&
&,ITER,ILHS,IRHS,A,B,C,D,E,F,G,H)
!
IF(TD.GE.TEND) GO TO 27
GO TO 28
!
! WRITE FINAL RESULTS
!
27 WRITE(6,96)N,LHORZ,LVERT,W,SIGMA,THETA,DT,DH,ERRSOR,ITMAX&
&,DNRM,ITER,TEND
WRITE(6,202)TD
!
95 FORMAT(7(TR1,E10.3))
96 FORMAT(//,' NO. OF INTERVALS PER LENGTH=',I3,/,TR1,'NO. OF&
&LENGTHS ON THE HORIZONTAL=',I4,/,TR1,'NO. OF LENGTHS ON THE&
&VERTICAL=',I4,/,TR1,'CONVERGE FACTOR=',F4.2,/,TR1,'SIGMA=',e11.4&
&,/,TR1,'THETA=',F4.2,/,TR1,'TIME INTERVAL=',E15.8,/,TR1,'STEP-SIZE IN&
&X AND Y CORDINATES=',E15.8,/,TR1,'ERRSOR=',E15.8,/,TR1,'ITMAX=',I4&
&,/,TR1,'MEASURE OF LARGEST ERROR BETWEEN LAST TWO ITERATES'&
&,E15.8,/,TR1,'NUMBER OF ITERATIONS DONE',I4,/,TR1,'FINAL TIME',e11.4&
&,///,TR1,'FINAL SOLUTIONS ARE STORED IN FILE CONCEN.ASCII ....',//)
!
202 FORMAT(TR1,'TIME=',E11.4,//)
203 FORMAT(TR1,' I-th TIME STEP=',I4,' TIME=',E11.4,' WT% @ 0.5mm =',E11.4)
204 FORMAT(TR1,I4,TR1,E11.4,TR1,E11.4)
STOP
CONTAINS
SUBROUTINE SLEEP(N,LHORZ,LVERT,W,U,ERRSOR,ITMAX,DNRM&
&,ITER,ILHS,IRHS,A,B,C,D,E,F,G,H)
IMPLICIT NONE
REAL, DIMENSION(:), INTENT(IN OUT) :: A, B, C, D, E, F, G, H, U
INTEGER, DIMENSION(:), INTENT(IN) :: ILHS, IRHS
INTEGER, INTENT(IN OUT) :: N, LHORZ, LVERT, ITER, ITMAX
INTEGER :: NM, IC, IT, J, K, IB, ID
REAL, INTENT(IN OUT) :: W, ERRSOR, DNRM
REAL :: T
NM=LVERT*N-1
IC=LHORZ*N+1
IT=LHORZ*N+3
DO ITER=1,ITMAX
DO J=1,NM
IB=ILHS(J)
ID=IRHS(J)
DO K=IB,ID
F(K)=G(K)*U(K)+H(K)*(U(K-IC)+U(K-1)+U(K+1)+U(K+IC))
T=U(K)+(W/C(K))*(F(K)-(A(K)*U(K-IC)+B(K)&
&*U(K-1)+C(K)*U(K)+D(K)*U(K+1)+E(K)*U(K+IC)))
IF(K.EQ.IT) GO TO 1
GO TO 2
1 DNRM=ABS(T-U(K))
GO TO 3
2 DNRM=MAX(DNRM,ABS(T-U(K)))
3 U(K)=T
END DO
END DO
IF(DNRM.LE.ERRSOR) GO TO 5
END DO
WRITE(6,6)
5 RETURN
6 FORMAT(' DID NOT CONVERGE IN SLEEP')
END SUBROUTINE SLEEP
SUBROUTINE ESTAB(NN,LHORZ,LVERT,SIGMA,THETA,DT,DH&
&,A,B,C,D,E,F,G,H)
IMPLICIT NONE
REAL, DIMENSION(:), INTENT(IN OUT) :: A, B, C, D, E, F, G, H
INTEGER, INTENT(IN) :: NN, LHORZ,LVERT
INTEGER :: J, NM, IC, M
REAL, INTENT(IN) :: SIGMA, THETA, DT, DH
REAL :: V, X, Y, Z
V=SIGMA*THETA*DT/DH**2
X=1+4*V
Y=SIGMA*(1-THETA)*DT/DH**2
Z=1-4*Y
NM=LVERT*NN-1
IC=LHORZ*NN+1
M=IC*(NM+2)
DO J=1,M
A(J)=-V
B(J)=-V
C(J)=X
D(J)=-V
E(J)=-V
F(J)=0.0
G(J)=Z
H(J)=Y
END DO
RETURN
END SUBROUTINE ESTAB
END PROGRAM DIFF2