#include<stdio.h>
#include<math.h>

#define dt   0.001                    /* 時間の刻み幅        */
#define dr   0.01                     /* 空間の刻み幅        */
#define imx  26                       /* x方向の最大格子点数 */
#define jmx  26                       /* y方向の最大格子点数 */
#define tmx  5000                     /* 最大 TimeStep 数    */

#define nu   1E-3                     /* 動粘性率            */
#define e    1E-5                     /* 収束判定条件        */
#define U    4.04E-4                  /* 上辺壁速度設定      */
#define buf  20                       /* 出力ファイル名buf   */
#define step 100                      /* ファイル出力間隔    */

#define imm (imx-1)
#define jmm (jmx-1)
#define imp (imx+1)
#define jmp (jmx+1)
#define hdr (dr /2)

static double pwr[imx+2][jmx+2];
static double npw[imx+2][jmx+2];
static double trs[imx+2][jmx+2];

static double  vx[imx+2][jmx+2];
static double  vy[imx+2][jmx+2];
static double  nx[imx+2][jmx+2];
static double  ny[imx+2][jmx+2];
static double hvx[imx+2][jmm+2];
static double hvy[imm+2][jmx+2];
static double mvx[imm+2][jmx+2];
static double mvy[imm+2][jmx+2];

static double  DD[imx+2][jmx+2];
static double  amp,XX,YY;
static char   fname[buf];
static int    i,j;
static int    lap = 0;

void init_value(void);

void solv_power(void);
void solv_half(void);
void solv_velocity(void);
void chan_cases(void);

void diag_bunch(void);
void diag_velocity(void);
void diag_power(void);

int main()
{
  init_value();
  diag_bunch();

  for(lap=1;lap<=tmx;lap++)
    {
      solv_half();

      solv_power();
      solv_velocity();

      chan_cases();
      diag_bunch();

      printf("%6.2f %% \n",(double)(lap/(double)tmx*100));
    }

  return 0;
}

void init_value(void)
{
  /* All parameter init */
  for(i=0;i<=imx+1;i++)
    for(j=0;j<=jmx+1;j++)
      pwr[i][j] = npw[i][j] = 0.0;

  for(i=0;i<=imx+1;i++)
    for(j=0;j<=jmx+1;j++)
      nx[i][j] = vx[i][j] = ny[i][j] = vy[i][j] = 0.0;

  for(i=0;i<=imx+1;i++)
    for(j=0;j<=jmx+1;j++)
      DD[i][j] = 0.0;
}

void solv_half(void)
{
  /* hvx[i][j] = vx(i+1/2,j+1/2); hvy[i][j] = vy(i+1/2,j+1/2) */
  for(i=0;i<=imx+1;i++)
    for(j=0;j<jmx;j++)
      hvx[i  ][j] = (vx[i][j+1] + vx[i][j])/2;

  for(i=0;i<imx;i++)
    for(j=0;j<=jmx+1;j++)
      hvy[i][j  ] = (vy[i+1][j] + vy[i][j])/2;

  /* mvx[i][j] = vx(i,j) */
  for(i=0;i<=imx;i++)
    for(j=0;j<=jmx;j++)
      mvx[i+1][j] = (vx[i+1][j] + vx[i][j])/2;

  /* mvy[i][j] = vy(i,j) */
  for(i=0;i<=imx;i++)
    for(j=0;j<=jmx;j++)
      mvy[i][j+1] = (vy[i][j+1] + vy[i][j])/2;
}

void solv_velocity(void)
{
  double vxx,vxy,vyy;
  double v2x,v2y,pxy;

  /* solve vx & vy */
  for(i=1;i<=imx;i++)
    for(j=1;j<jmx;j++)
      {
        vxx = mvx[i][j]*mvx[i  ][j  ]-mvx[i-1][j]*mvx[i-1][j];
        vxy = hvx[i][j]*hvy[i-1][j+1]-hvx[i][j-1]*hvy[i-1][j];

        v2x = (vx[i+1][j]-2*vx[i][j]+vx[i-1][j])/(dr*dr);
        v2y = (vx[i][j+1]-2*vx[i][j]+vx[i][j-1])/(dr*dr);

        pxy = (pwr[i][j]-pwr[i-1][j])/dr;

        XX  = vxx/dr + vxy/dr;
        nx[i][j] = vx[i][j] + (XX-pxy+(v2x+v2y)*nu)*dt;
      }

  for(i=1;i<imx;i++)
    for(j=1;j<=jmx;j++)
      {
        vyy = mvy[i][j]*mvy[i  ][j  ]-mvy[i][j-1]*mvy[i][j-1];
        vxy = hvx[i][j]*hvy[i+1][j-1]-hvx[i-1][j]*hvy[i][j-1];

        v2x = (vy[i+1][j]-2*vy[i][j]+vy[i-1][j])/(dr*dr);
        v2y = (vy[i][j+1]-2*vy[i][j]+vy[i][j-1])/(dr*dr);

        pxy = (pwr[i][j]-pwr[i][j-1])/dr;

        YY  = vyy/dr + vxy/dr;
        ny[i][j] = vy[i][j] + (YY-pxy+(v2x+v2y)*nu)*dt;
      }

  /* Boundary Conditions */
  /* Free Wall --side */
  for(j=1;j<jmx;j++)
    {
      nx[  1][j] =  nx[imx][j] = 0.0;

      nx[  0][j] = -nx[  2][j];
      nx[imp][j] = -nx[imm][j];

      ny[  0][j] =  ny[  1][j];
      ny[imx][j] =  ny[imm][j];
    }

  for(i=1;i<imx;i++)
    {
      ny[i][  1] =  ny[i][jmx] = 0.0;

      ny[i][  0] = -ny[i][  2];
      ny[i][jmp] = -ny[i][jmm];

      nx[i][  0] =  nx[i][  1];
 
      nx[i][jmx] =  nx[i][jmm] + amp;
    }
}

void solv_power(void)
{
  double v2x,v2y,vxy;
  double err,sum;

  for(i=0;i<=imx+1;i++)
   for(j=0;j<=jmx+1;j++)
     npw[i][j] = pwr[i][j];

  for(i=1;i<imx;i++)
    for(j=1;j<jmx;j++)
      DD[i][j] = (vx[i+1][j]-vx[i][j]+vy[i][j+1]-vy[i][j])/dr;

  for(i=1;i<imx;i++)
    for(j=1;j<jmx;j++)
      {
        v2x = pow(mvx[i+1][j],2)-2*pow(mvx[i][j],2)+pow(mvx[i-1][j],2);
        v2y = pow(mvy[i][j+1],2)-2*pow(mvy[i][j],2)+pow(mvy[i][j-1],2);

        vxy = hvx[i+1][j]*hvy[i][j+1] - hvx[i+1][j-1]*hvy[i][j]
          - hvx[i][j]*hvy[i-1][j+1] + hvx[i][j-1]*hvy[i-1][j];

        trs[i][j] = (v2x + 2*vxy + v2y) - DD[i][j]*dr*dr/dt
          -(DD[i+1][j]+DD[i-1][j]-4*DD[i][j]+DD[i][j+1]+DD[i][j-1])*nu;
      }

  /* Gauss - Seidel Method */
  for(;;)
    {
      err = 0.0;
      for(i=1;i<imx;i++)
        {
          for(j=1;j<jmx;j++)
            {
              sum = npw[i][j];

              npw[i][j] = ((npw[i][j-1]+npw[i][j+1]+npw[i-1][j]+npw[i+1][j])
                           +trs[i][j  ])/4;

              err += fabs(npw[i][j]-sum);
            }
        }
      /* Free Wall --side */
      for(j=1;j<jmx;j++)
        {
          npw[  0][j] = npw[  1][j];
          npw[imx][j] = npw[imm][j];
        }
      for(i=1;i<imx;i++)
        {
          npw[i][  0] = npw[i][  1];
          npw[i][jmx] = npw[i][jmm];
        }
      if(err <= e) break;
    }

  /* Moving Wall conditions */
  amp = U;
}

void chan_cases(void)
{
  for(i=0;i<=imx+1;i++)
    for(j=0;j<=jmx+1;j++)
      pwr[i][j] = npw[i][j];

  for(i=0;i<=imx+1;i++)
    for(j=0;j<=jmx+1;j++)
      vx[i][j] = nx[i][j];

  for(i=0;i<=imx+1;i++)
    for(j=0;j<=jmx+1;j++)
      vy[i][j] = ny[i][j];
}

void diag_bunch(void)
{
  if(lap % step == 0)
    {
      diag_power();
      diag_velocity();
    }
}

void diag_power(void)
{
  FILE *fp;

  sprintf(fname,"pwr%05d.dat",lap);
  fp = fopen(fname,"w");

  for(i=1;i<imx;i++)
    {
      for(j=1;j<jmx;j++)
        fprintf(fp,"%e %e %e\n",(i-1)*dr,(j-1)*dr,pwr[i][j]);
      fprintf(fp,"\n");
    }

  fclose(fp);
}

void diag_velocity(void)
{
  double avx,avy;
  FILE *fp;

  sprintf(fname,"vxy%05d.dat",lap);
  fp = fopen(fname,"w");

  for(i=1;i<imx;i++)
    {
      for(j=1;j<jmx;j++)
        {
          avx = (vx[i+1][j]+vx[i+1][j+1])/2*50;
          avy = (vy[i][j+1]+vy[i+1][j+1])/2*50;
          fprintf(fp,"%e %e %e %e\n",(i-1)*dr,(j-1)*dr,avx,avy);
        }
      fprintf(fp,"\n");
    }

  fclose(fp);
}