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)) ; }
  }