//=================================================================C
//  2-Dimensional Cavity Flow by velocity and pressure             C
//       with FDM and SOR                                          C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,   2007/06/04                     C
//        ( Tokyo Polytechnic University )                         C
//    後 保範（東京工芸大学）                                      C
//=================================================================C
#include <stdio.h>
#include <math.h>
#include <time.h>
// Global Define
#define NDX  301
#define NDY  301
double P[NDX][NDY], B[NDX][NDY] ;
double U[NDX][NDY], V[NDX][NDY] ;
double UW[NDX][NDY], VW[NDX][NDY] ;
FILE   *FT1 ;
// Function Definition
int PSOR(int NX, int NY, double DT, double H, double ALP, double EPS) ;
void OUTUVP(int NX, int NY) ;
// Main Program
int main()
{ int     i, j, k, NX, NY, NT, MT, ID, T1, T2 ;
  int     ITER, NTER, NI1, NI2, IER ;
  double  Re, DT, T, H, DH, H2, D2, D5, EPS, ALP, CPU, OVER ;
  char    FLOW2[13]={'F','L','O','W','2','-','?','?','.','t','x','t',} ;
  char    NTABL[10]={'0','1','2','3','4','5','6','7','8','9'} ;
//  Initial Data
  ID  = 0 ;
  IER = 0 ;
  printf("Type In  NX,Re \n") ;
  scanf("%d %lf",&NX,&Re) ; 
  printf("DT,NT(Total NO.),MT(Output Interval) \n") ;
  scanf("%lf %d %d",&DT,&NT,&MT) ;
  if(NX >= NDX) { NX = NDX - 1 ; }
  NY = NX ;
  printf("# Cavity Flow Numerical Analysis \n") ;
  printf("NX=%d NY=%d Re=%f DT=%f  NT,MT=%d %d \n",NX,NY,Re,DT,NT,MT) ;
//  Initial Constant
  T  = 0.0 ;
  DH = NX ;
  D5 = DH/2.0 ;
  H  = 1.0/DH ;
  H2 = H*H ;
  D2 = DH*DH ;
//  Set Initial Condition
  for (i=0; i<=NX; i++) {
    for (j=0; j<=NY; j++) {
      P[i][j] = 0.0 ;
      U[i][j] = 0.0 ;
      V[i][j] = 0.0 ;  }
    }
//   U=1.0 on L1
  for (i=0; i<=NX; i++) {
    U[i][NY] = 1.0 ;  }
//  SOR Parameter (ALP) 
  ALP = 1.0 + log(NX*1.0)/log(NDX*1.2) ;
  EPS = 1.0e-4 ;
//  Main Loop
  NTER = 0 ;
  T1 = clock() ;
  for (k=1; k<=NT; k++) {
    T = T + DT ;
//   Compute U,V 
    OVER = 0.0 ;
    for (i=1; i<NX; i++) {
      for (j=1; j<NY; j++) {
//        U[i][j] = U[i][j] + DT*((U[i-1][j]*U[i-1][j]-U[i+1][j]*U[i+1][j]
          UW[i][j] = U[i][j] + DT*((U[i-1][j]*U[i-1][j]-U[i+1][j]*U[i+1][j]
                + U[i][j-1]*V[i][j-1] - U[i][j+1]*V[i][j+1]
                + P[i-1][j] - P[i+1][j])*D5 - (4.0*U[i][j] - U[i][j-1]
                - U[i-1][j] - U[i+1][j] - U[i][j+1])*D2/Re ) ;
//        V[i][j] = V[i][j] + DT*((V[i][j-1]*V[i][j-1]-V[i][j+1]*V[i][j+1]
          VW[i][j] = V[i][j] + DT*((V[i][j-1]*V[i][j-1]-V[i][j+1]*V[i][j+1]
                + U[i-1][j]*V[i-1][j] - U[i+1][j]*V[i+1][j]
                + P[i][j-1] - P[i][j+1])*D5 - (4.0*V[i][j] - V[i][j-1]
                - V[i-1][j] - V[i+1][j] - V[i][j+1])*D2/Re ) ;
        OVER = OVER + fabs(U[i][j]) + fabs(V[i][j]) ;  }
       } 
//   Copy Back
    for (i=1; i<NX; i++) {
      for (j=1; j<NY; j++) {
        U[i][j] = UW[i][j] ;
        V[i][j] = VW[i][j] ; }
	  }
//   Check over-flow computation
    if(OVER >= 1.0e10) {
      printf("** Stop for over-flow computation ** \n") ;
      printf(" You have to give smaller DT \n") ;
      IER = 1 ;
      goto M100 ;  }
//   Compute P
    ITER = PSOR(NX, NY, DT, H, ALP, EPS) ;
    NTER = NTER + ITER ;
//   Output U,V,Phi
    if( k%MT == 0) {
      NI1      = ID/10 ;
      NI2      = ID - NI1*10 ;
      FLOW2[6] = NTABL[NI1] ;
      FLOW2[7] = NTABL[NI2] ;
      FT1 = fopen(FLOW2,"w") ;
      fprintf(FT1,"# Cavity Flow Numerical Analysis \n") ;
      fprintf(FT1,"NX=%d NY=%d Re=%f DT=%f NT,MT=%d %d \n",NX,NY,Re,DT,NT,MT) ;
      fprintf(FT1,"#Time=%7.2f,  Step=%5d  SOR-loop=%6d \n",T,k,NTER) ;
      printf("#Time=%7.2f,  Step=%5d SOR-loop=%6d \n",T,k,NTER) ;
      NTER = 0 ;
      OUTUVP(NX, NY) ;
      fclose(FT1) ;
      if(ID < 99) { ID = ID + 1 ; }
     }
   }
  T2 = clock() ;
  CPU  = (double)(T2 - T1)/CLOCKS_PER_SEC ;
  printf(" NX=%4d,  NT=%5d,  Time(s)=%9.2f \n",NX,NT,CPU) ;
  M100:  return (IER) ;
  } 
//=================================================================C
int PSOR(int NX, int NY, double DT, double H, double ALP, double EPS)
//=================================================================C
//  Solve Ax=b by SOR with 2 dimensional FDM                       C
//    Get P with U,V ( Acceleration factor )                       C
//-----------------------------------------------------------------C
//    U,V,P,B ;  Global Array                                      C
//    NX      I*4, In,  Grid Numbers on X-axis                     C
//    NY      I*4, In,  Grid Numbers on Y-axis                     C
//    DT      R*8, In,  Delta T                                    C
//    H       R*8, In,  H=1.0/NX                                   C
//    ALP     R*8, In,  SOR Acceleration factor                    C
//    EPS     R*8, In,  if ||r||/||b|| <= EPS --> return           C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,   2007/06/04                     C
//        ( Tokyo Polytechnic University )                         C
//=================================================================C
{ int     i, j, k, IP, IM, JP, JM, ITER ;
  double  BN, RN, ERR, R ;
  double  UP0, UM0, U0P, U0M, VP0, VM0, V0P, V0M ;
//  Set B and BN (B norm)
  BN = 0.0 ;
  for (i=0; i<=NX; i++) {
    for (j=0; j<=NY; j++) {
//   Boundary on L1
      if(j != NY) { U0P=U[i][j+1] ;  V0P=V[i][j+1] ;  }
        else { U0P=1.0 ;  V0P=0.0 ;  } 
//   Boundary on L2
      if(i != 0)  { UM0=U[i-1][j] ;  VM0=V[i-1][j] ;  }
        else { UM0=U[1][j] ;  VM0=-V[1][j] ;  }
//   Boundary on L3
      if(j != 0)  { U0M=U[i][j-1] ;  V0M=V[i][j-1] ;  }
        else { U0M=-U[i][1] ;  V0M=V[i][1] ;  }
//   Boundary on L4
      if(i != NX) { UP0=U[i+1][j] ;  VP0=V[i+1][j] ;  }
        else { UP0=U[NX-1][j] ;  VP0=-V[NX-1][j] ;  }
//   Computation
      B[i][j] = ((UP0-UM0)*(UP0-UM0) + (V0P-V0M)*(V0P-V0M))/4.0
              + (VP0 - VM0)*(U0P - U0M)/2.0
              - (UP0 - UM0 + V0P - V0M)*H/(2.0*DT) ;
      BN = BN + B[i][j]*B[i][j] ;  }
     }
//  Main Loop
  for (k=1; k<=NX*NY; k++) {
    RN = 0.0 ;
    for (i=0; i<=NX; i++) {
      IP = i + 1 ;
      IM = i - 1 ;
//   Boundary on L1,L3
      if(i == NX) { IP = NX - 1 ; }
      if(i == 0)  { IM = 1 ;  }
      for (j=0; j<=NY; j++) {
        JP = j + 1 ;
        JM = j - 1 ;
//   Boundary on L2,L4
        if(j == NY) { JP = NY - 1 ; }
        if(j == 0)  { JM = 1 ; }
//   Computation
         R = (B[i][j]+P[i][JM]+P[IM][j]+P[IP][j]+P[i][JP])/4.0-P[i][j] ;
         P[i][j] = P[i][j] + ALP*R ;
//   Reset on P1
         if(i==0 & j==0) { P[i][j]=0.0 ;  R=0.0 ;  }
         RN = RN + R*R*16.0 ;
     }  }
//   if(ERR <= EPS) return
    ERR = sqrt(RN/BN) ;
    if(ERR <= EPS) goto M100 ;
    }
  M100:  ITER = k ;
  return (ITER) ;
  }
//=================================================================C
void OUTUVP(int NX, int NY)
//=================================================================C
//  Output U,V,P                                                   C
//-----------------------------------------------------------------C
//    P,U,V ; Global Array                                         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,   2007/06/04                     C
//        ( Tokyo Polytechnic University )                         C
//=================================================================C
{ int     i, j ;
  double  X, Y, DX, DY ;
//   Output
  DX = 1.0/NX ;
  DY = 1.0/NY ;
  for (j=0; j<=NY; j++) {
    Y = j*DY ;
    for (i=0; i<=NX; i++) {
      X = i*DX ;
      fprintf(FT1,"%9.5f %9.5f %9.5f %9.5f %9.5f \n",
                  X,Y,P[i][j],U[i][j],V[i][j]) ;  }
    }
  }