3次元FDM用SOR法(FORTRAN)テストプログラム
2003/08/24 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
3次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してSOR法で
反復解xを求めるプログラムのテストプログラム。
加速係数(ω)を与える必要がある。収束速度はωにより大きく左右される。
2. プログラム
C=================================================================C
C Test Program of SOR for Dimensional FDM C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/08/17 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
PARAMETER(NDX=50, NDY=50, NDZ=50, NDXY=(NDX+1)*(NDY+1)*(NDZ+1))
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(7*NDXY), B(NDXY), X(NDXY)
C Initial Data
OPEN(UNIT=1,FILE='F-SOR3D.txt')
CALL XCLOCK(T1,3)
EPS = 1.0e-6
WRITE(6,*) 'Type In NX,NY,NZ,Omega'
READ(5,*) NX,NY,NZ,Omega
if(NX.gt.NDX) NX = NDX
if(NY.gt.NDY) NY = NDY
if(NZ.gt.NDZ) NZ = NDZ
WRITE(6,100) NX,NY,NZ,Omega
WRITE(1,100) NX,NY,NZ,Omega
C Set A,B,X
CALL SET3D(A,B,X,NX,NY,NZ)
C Solve Ax=b by SOR
ITER = 5*NX*NY*NZ
CALL XCLOCK(T1,5)
CALL SOR3D(A,B,NX,NY,NZ,X,Omega,EPS,ITER,ERR,IERR)
F1 = 2
F2 = 18
CALL XCLOCK(T2,5)
C Write
CPU = T2 - T1
FLOP = (NX-1)*(NY-1)*(NZ-1)*(F1+F2*ITER)*1.0D-6
FLOPS = 0.0
if(CPU.ne.0.0) FLOPS = FLOP/CPU
WRITE(6,200) ITER,ERR,CPU,FLOPS
WRITE(1,200) ITER,ERR,CPU,FLOPS
C Output center-line
CALL OUT3D(X,NX,NY,NZ)
C
STOP
100 FORMAT(1H ,'SOR Method for 2-Dimensional FDM, NX,NY,NZ,Omega=',
13I3,F6.3,/,' Loop Error Time(s) MFLOPS')
200 FORMAT(1H ,I6,E10.2,F8.2,F8.1)
END
C=================================================================C
SUBROUTINE SET3D(A,B,X,NX,NY,NZ)
C=================================================================C
C Set A and B by 3 Dimensional FDM C
C -div(K.grad(X)) = f, K=1.0, f=10,0 and 1x1 square C
C-----------------------------------------------------------------C
C A(0:NX,0:NY,0:NZ,7) R*8, Out, A Coefficient Matrix C
C B(0:NX,0:NY,0:NZ) R*8, Out, A Right-hand Vector(b) C
C X(0:NX,0:NY,0:NZ) R*8, Out, Initial X, X=0 C
C NX I*4, In, Grid Numbers on X-axis C
C NY I*4, In, Grid Numbers on Y-axis C
C NZ I*4, In, Grid Numbers on Z-axis C
C ID I*4, In, ID=0 ; NO output Un-Symmetric C
C =1 ; Output Un-symmetric C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/08/17 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(0:NX,0:NY,0:NZ,7), B(0:NX,0:NY,0:NZ)
DIMENSION X(0:NX,0:NY,0:NZ)
C Clear A,B,X
do k=0,NZ
do j=0,NY
do i=0,NX
A(i,j,k,1) = 0.0
A(i,j,k,2) = 0.0
A(i,j,k,3) = 0.0
A(i,j,k,4) = 0.0
A(i,j,k,5) = 0.0
A(i,j,k,6) = 0.0
A(i,j,k,7) = 0.0
B(i,j,k) = 0.0
X(i,j,k) = 0.0
end do
end do
end do
C Initial Data
DX = NX
DY = NY
DZ = NZ
HX = 1.0/DX
HY = 1.0/DY
HZ = 1.0/DZ
F = 100.0
C Set A,B
do k=1,NZ-1
do j=1,NY-1
do i=1,NX-1
A(i,j,k,1) = -HX*HY*DZ
if(k.eq.1) A(i,j,k,1) = 0.0
A(i,j,k,2) = -HZ*HX*DY
if(j.eq.1) A(i,j,k,2) = 0.0
A(i,j,k,3) = -HY*HZ*DX
if(i.eq.1) A(i,j,k,3) = 0.0
A(i,j,k,4) = 2.0*(HY*HZ*DX+HZ*HX*DY+HX*HY*DZ)
A(i,j,k,5) = -HY*HZ*DX
if(i.eq.NX-1) A(i,j,k,5) = 0.0
A(i,j,k,6) = -HZ*HX*DY
if(j.eq.NY-1) A(i,j,k,6) = 0.0
A(i,j,k,7) = -HX*HY*DZ
if(k.eq.NZ-1) A(i,j,k,7) = 0.0
B(i,j,k) = HX*HY*HZ*F
end do
end do
end do
C
RETURN
END
C=================================================================C
SUBROUTINE OUT3D(X,NX,NY,NZ)
C-----------------------------------------------------------------C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION X(0:NX,0:NY,0:NZ)
C Output Y-Z center Line
WRITE(1,100)
DX = NX
do i=0,NX
XV = i/DX
WRITE(1,110) XV,X(i,NY/2,NZ/2)
end do
C
100 FORMAT(1H ,/'Output on Y,Z-center line'/' X-axis Solution(X)')
110 FORMAT(1H ,2F10.6)
RETURN
END
C=================================================================C
SUBROUTINE XCLOCK(CPU,ID)
C=================================================================C
C CPU time Subroutine C
C CPU R*8 Out, CPU Time C
C ID I*4 In, Dummy ( Same Hitachi FORTRAN XCLOCK ) C
C-----------------------------------------------------------------C
C Written by Y. Ushiro ( Waseda and Hitachi ) 2002.10.29 C
C=================================================================C
REAL*8 CPU
INTEGER*2 I1, I2, I3, I4
C
IF(ID.GE.1) THEN
CALL GETTIM(I1,I2,I3,I4)
CPU = ( I1*60.0 + I2 )*60.0 + I3 + I4*0.01
END IF
C
RETURN
END