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