3次元FDM用SOR法(C)テストプログラム
2003/08/24 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
3次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してSOR法で
反復解xを求めるプログラムのテストプログラム。
加速係数(ω)を与える必要がある。収束速度はωにより大きく左右される。
2. プログラム
//=================================================================C
// Test Program of SOR for 3-Dimensional FDM C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/08/21 C
// ( Hitachi Ltd. and Waseda University ) C
//=================================================================C
#include <stdio.h>
#include <math.h>
#include <time.h>
// Global Define
#define NDX 51
#define NDY 51
#define NDZ 51
#define A(i,j,k,l) A[l-1][k][j][i]
#define B(i,j,k) B[k][j][i]
#define X(i,j,k) X[k][j][i]
double A[7][NDZ][NDY][NDX], B[NDZ][NDY][NDX], X[NDZ][NDY][NDX] ;
FILE *FT1 ;
// Function Definition
void SET3D(int NX, int NY, int NZ) ;
void OUT3D(int NX, int NY, int NZ) ;
int SOR3D(int NX, int NY, int NZ, double Omega, double EPS,
int *ITER,double *ERR) ;
// Main Program
void main ()
{ int NX, NY, NZ, ITER, IERR, T1, T2 ;
double EPS, ERR, CPU, F1, F2, FLOP, FLOPS, Omega ;
// Open File
FT1 = fopen("C-SOR3D.txt","w") ;
EPS = 1.0e-6 ;
// Size Input
printf("Type In NX,NY,NZ,Omega \n") ;
scanf("%d %d %d %lf",&NX,&NY,&NZ,&Omega) ;
if(NX > NDX-1) NX = NDX-1 ;
if(NY > NDY-1) NY = NDY-1 ;
if(NZ > NDZ-1) NZ = NDZ-1 ;
fprintf(FT1,"SOR Method for 3-Dimensional FDM") ;
fprintf(FT1," NX=%d, NY=%d NZ=%d Omega=%f \n",NX,NY,NZ,Omega) ;
fprintf(FT1," Loop Error Time(s) MFLOPS \n") ;
// Set A,B,X
SET3D(NX, NY, NZ) ;
// Solve Ax=b by CG
ITER = 5*NX*NY*NZ ;
T1 = clock() ;
IERR = SOR3D(NX, NY, NZ, Omega, EPS, &ITER, &ERR) ;
F1 = 2 ;
F2 = 18 ;
T2 = clock() ;
// Write
CPU = (double)(T2 - T1)/CLOCKS_PER_SEC ;
FLOP = (NX-1)*(NY-1)*(NZ-1)*(F1+F2*ITER)*1.0e-6 ;
FLOPS = 0.0 ;
if(CPU > 0.0) FLOPS = FLOP/CPU ;
fprintf(FT1,"%6d %10.2e %8.1f %8.1f \n",ITER,ERR,CPU,FLOPS) ;
printf("%6d %10.2e %8.1f %8.1f \n",ITER,ERR,CPU,FLOPS) ;
// Output center-line
OUT3D(NX, NY, NZ) ;
}
//=================================================================C
void SET3D(int NX, int NY, int NZ)
//=================================================================C
// Set A and B by 3 Dimensional FDM with symmetric C
// -div(K.grad(X)) = f, K=1.0, f=10,0 and 1x1 square C
//-----------------------------------------------------------------C
// NX I*4, In, Grid Numbers on X-axis C
// NY I*4, In, Grid Numbers on Y-axis C
// NZ I*4, In, Grid Numbers on Z-axis C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/08/21 C
// ( Hitachi Ltd. and Waseda University ) C
//=================================================================C
{ int i, j, k ;
double DX, DY, DZ, HX, HY, HZ, F ;
// Clear A,B,X
for (k=0; k<=NY; k++) {
for (j=0; j<=NY; j++) {
for (i=0; i<=NX; i++) {
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 ; }
} }
// Initial Data
DX = NX ;
DY = NY ;
DZ = NZ ;
HX = 1.0/DX ;
HY = 1.0/DY ;
HZ = 1.0/DZ ;
F = 100.0 ;
// Set A,B
for (k=1; k<=NZ-1; k++) {
for (j=1; j<=NY-1; j++) {
for (i=1; i<=NX-1; i++) {
A(i,j,k,1) = -HX*HY*DZ ;
if(k == 1) A(i,j,k,1) = 0.0 ;
A(i,j,k,2) = -HZ*HX*DY ;
if(j == 1) A(i,j,k,2) = 0.0 ;
A(i,j,k,3) = -HY*HZ*DX ;
if(i == 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 == NX-1) A(i,j,k,5) = 0.0 ;
A(i,j,k,6) = -HZ*HX*DY ;
if(j == NY-1) A(i,j,k,6) = 0.0 ;
A(i,j,k,7) = -HX*HY*DZ ;
if(k == NZ-1) A(i,j,k,7) = 0.0 ;
B(i,j,k) = HX*HY*HZ*F ; }
} }
}
//=================================================================C
void OUT3D(int NX, int NY, int NZ)
//-----------------------------------------------------------------C
{ int i ;
double DX, XV ;
// Output Y
fprintf(FT1,"Output on YZ-center line \n"," X-axis Solution(X) \n") ;
DX = NX ;
for (i=0; i<=NX; i++) {
XV = i/DX ;
fprintf(FT1,"%10.6f %10.6f \n",XV,X(i,NY/2,NZ/2)) ; }
}