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