// cavity flow
// The source code is written in C programming language.
// Copyright @2026 Takeshi Seta All Rights Reserved.
//
// 5  1  6
//    |
// 3--0--4
//    |
// 7  2  8
//
// rho  : density
// u    : horizontal velocity component (time = n)
// v    : vertical velocity component (time = n)
// un   : horizontal velocity component (time = n -1)
// vn   : vertical velocity component (time = n - 1)
// ut   : horizontal velocity component of the top wall
// vt   : vertical velocity component of the top wall
// nx   : number of grid points (x-axis)
// ny   : number of grid points (y-axis)
// f    : distribution function
// f0   : equilibrium distribution function
// cx   : discrete velocity (x-axis)
// cy   : discrete velocity (y-axis)
// tau  : relaxation time
// mu   : kinematic viscosity
// re   : Reynolds number
// omega: vorticity

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

# define DIM 103 

int main(void)
{
  FILE   *fp;
  int    nx = 101, ny = 101, time = 0, loop1, loop2;
  int    i, j, k, in, jn, ip, im, jp, jm;
  double rho[DIM][DIM], u[DIM][DIM], v[DIM][DIM], un[DIM][DIM], vn[DIM][DIM], omega[DIM][DIM];
  double f[9][DIM][DIM], f0[9][DIM][DIM], ftmp[9][DIM][DIM], cx[9], cy[9];
  double ut = 0.050, vt = 0.0, umax, umin, tmp, u2, mu, norm, err, tau, re = 1000;
  char   a[DIM][DIM], filename[256];
  FILE *fp_pvd;
  fp_pvd = fopen("lbmcavi.pvd", "w");

  fprintf(fp_pvd, "<VTKFile type=\"Collection\" version=\"0.1\">\n");
  fprintf(fp_pvd, "  <Collection>\n");

  // initial condition
  mu = ut*(double)(nx - 1)/re;
  tau = 3.0*mu + 0.5;

  for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){
    u[i][j] = 0.0; v[i][j] = 0.0; un[i][j] = 0.0; vn[i][j] = 0.0;
    rho[i][j] = 1.0;
  } }

  cx[0] = 0.0; cy[0] = 0.0;
  cx[1] = 0.0; cy[1] = 1.0;
  cx[2] = 0.0; cy[2] =-1.0;
  cx[3] =-1.0; cy[3] = 0.0;
  cx[4] = 1.0; cy[4] = 0.0;
  cx[5] =-1.0; cy[5] = 1.0;
  cx[6] = 1.0; cy[6] = 1.0;
  cx[7] =-1.0; cy[7] =-1.0;
  cx[8] = 1.0; cy[8] =-1.0;

  for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){
    u2 = u[i][j]*u[i][j] + v[i][j]*v[i][j];      
    f0[0][i][j] = rho[i][j]*(1.0 -3.0/2.0*u2)*4.0/9.0;
    for(k = 1; k <= 4; k++){
      tmp = cx[k]*u[i][j] + cy[k]*v[i][j];      
      f0[k][i][j] = rho[i][j]*(1.0 +3.0*tmp +9.0/2.0*tmp*tmp -3.0/2.0*u2)/9.0;
    }
    for(k = 5; k <= 8; k++){
      tmp = cx[k]*u[i][j] + cy[k]*v[i][j];      
      f0[k][i][j] = rho[i][j]*(1.0 +3.0*tmp +9.0/2.0*tmp*tmp -3.0/2.0*u2)/36.0;
    }
  } }

  for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){ for(k = 0; k <= 8; k++){
    f[k][i][j] = f0[k][i][j];
  } } }

// calculation start
  for(loop1 = 0; loop1 < 20; loop1++){
  for(loop2 = 0; loop2 < 5000; loop2++){
    time++;

    for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){
      un[i][j] = u[i][j]; vn[i][j] = v[i][j];
    } }

  // collision
    for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){
      u2 = u[i][j]*u[i][j] + v[i][j]*v[i][j];      
      f0[0][i][j] = rho[i][j]*(1.0 -3.0/2.0*u2)*4.0/9.0;
      for(k = 1; k <= 4; k++){
        tmp = cx[k]*u[i][j] + cy[k]*v[i][j];      
        f0[k][i][j] = rho[i][j]*(1.0 +3.0*tmp +9.0/2.0*tmp*tmp -3.0/2.0*u2)/9.0;
      }  
      for(k = 5; k <= 8; k++){
        tmp = cx[k]*u[i][j] + cy[k]*v[i][j];      
        f0[k][i][j] = rho[i][j]*(1.0 +3.0*tmp +9.0/2.0*tmp*tmp -3.0/2.0*u2)/36.0;
      }
    } }

    for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){ for(k = 0; k <= 8; k++){
      f[k][i][j] = f[k][i][j] - (f[k][i][j] - f0[k][i][j])/tau;
    } } }

    // propagation 
    for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){ for(k = 0; k <= 8; k++){
      ftmp[k][i][j] = f[k][i][j];
    } } }

    for(j = 0; j <= ny; j++){
      for(i = 1; i <= nx  ; i++){ f[3][i-1][j  ] = ftmp[3][i][j]; }
      for(i = 0; i <= nx-1; i++){ f[4][i+1][j  ] = ftmp[4][i][j]; }
    }

    for(j = 0; j <= ny-1; j++){
      for(i = 0; i <= nx  ; i++){ f[1][i  ][j+1] = ftmp[1][i][j]; }
      for(i = 1; i <= nx  ; i++){ f[5][i-1][j+1] = ftmp[5][i][j]; }
      for(i = 0; i <= nx-1; i++){ f[6][i+1][j+1] = ftmp[6][i][j]; }
    }
    for(j = 1; j <= ny; j++){
      for(i = 0; i <= nx  ; i++){ f[2][i  ][j-1] = ftmp[2][i][j]; }
      for(i = 1; i <= nx  ; i++){ f[7][i-1][j-1] = ftmp[7][i][j]; }
      for(i = 0; i <= nx-1; i++){ f[8][i+1][j-1] = ftmp[8][i][j]; }
    }

    // boundary condition
    // Half-way bounce back
    for(j = 0; j <= ny; j++){
      jm = j - 1; if(j ==  0){jm = ny;}
      jp = j + 1; if(j == ny){jp =  0;}
      f[4][   1][j] = f[3][0][j ];
      f[6][   1][j] = f[7][0][jm];
      f[8][   1][j] = f[5][0][jp];

      f[3][nx-1][j] = f[4][nx][j ];
      f[7][nx-1][j] = f[6][nx][jp];
      f[5][nx-1][j] = f[8][nx][jm];
    }

    for(i = 0; i <= nx; i++){
      im = i - 1; if(i ==  0){im = nx;}
      ip = i + 1; if(i == nx){ip =  0;}
      f[1][i][   1] = f[2][i ][ 0];
      f[6][i][   1] = f[7][im][ 0];
      f[5][i][   1] = f[8][ip][ 0];

      rho[i][ny] = (f[0][i][ny] + f[4][i][ny] + f[3][i][ny]
             + 2.0*(f[1][i][ny] + f[5][i][ny] + f[6][i][ny]))/(1.0 + vt);
      f[2][i][ny-1] = f[1][i ][ny] - rho[i][ny]*(cx[1]*ut + cy[1]*vt)*2.0/3.0;
      f[7][i][ny-1] = f[6][ip][ny] - rho[i][ny]*(cx[6]*ut + cy[6]*vt)/6.0;
      f[8][i][ny-1] = f[5][im][ny] - rho[i][ny]*(cx[5]*ut + cy[5]*vt)/6.0;
    }

    // corner (rho = 1.0)
    u2 = ut*ut + vt*vt;      
    tmp = cx[5]*ut + cy[5]*vt;      
    f[5][nx-1][ny-1] = (1.0 +3.0*tmp +9.0/2.0*tmp*tmp -3.0/2.0*u2)/36.0;
    tmp = cx[8]*ut + cy[8]*vt;      
    f[8][nx-1][ny-1] = (1.0 +3.0*tmp +9.0/2.0*tmp*tmp -3.0/2.0*u2)/36.0;

    tmp = cx[6]*ut + cy[6]*vt;      
    f[6][   1][ny-1] = (1.0 +3.0*tmp +9.0/2.0*tmp*tmp -3.0/2.0*u2)/36.0;
    tmp = cx[7]*ut + cy[7]*vt;      
    f[7][   1][ny-1] = (1.0 +3.0*tmp +9.0/2.0*tmp*tmp -3.0/2.0*u2)/36.0;

    f[6][nx-1][   1] = 1.0/36.0; f[7][nx-1][   1] = 1.0/36.0;
    f[5][   1][   1] = 1.0/36.0; f[8][   1][   1] = 1.0/36.0;

    // physics
    for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){
      rho[i][j] = f[0][i][j]; u[i][j] = 0; v[i][j] = 0;
      for( k = 1; k <= 8; k++){
        rho[i][j] = rho[i][j] + f[k][i][j];
          u[i][j] =   u[i][j] + f[k][i][j]*cx[k];
          v[i][j] =   v[i][j] + f[k][i][j]*cy[k];
      } 
      u[i][j] /= rho[i][j];
      v[i][j] /= rho[i][j];
    } }

    norm = 0.0;
    for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){
      tmp = sqrt(pow(u[i][j]-un[i][j],2) + pow(v[i][j]-vn[i][j],2));
      if(tmp > norm){ norm = tmp; }
    } }

  } //loop2
  
  tmp = 0.0;
  for(i = 1; i <= nx - 1; i++){ for(j = 1; j <= ny - 1; j++){
    tmp = tmp + rho[i][j];
  } }
  tmp = tmp/(nx - 1)/(ny - 1);

  err = 0.0;
  for(i = 1; i <= nx - 1; i++){ for(j = 1; j <= ny - 1; j++){
    err = err + pow((rho[i][j] - tmp),2);
  } }
  err = sqrt(err)/(nx - 1)/(ny - 1);
  err = err/tmp;

  printf("Time = %d, Norm = %8.6e, Error = %8.6e\n", time, norm, err);

  for(i = 1; i <= nx - 1; i++){ for(j = 1; j <= ny - 1; j++){
    a[i][j]='0';
  } }

  umax = -999.9; umin = 999.9;
  for(i = 1; i <= nx - 1; i++){ for(j = 1; j <= ny - 1; j++){
    tmp = sqrt(pow(u[i][j],2) + pow(v[i][j],2));
    if(tmp >= umax){umax = tmp;}
    if(tmp <= umin){umin = tmp;}
  } }

  for(i = 1; i <= nx - 1; i++){ for(j = 1; j <= ny - 1; j++){
    tmp = sqrt(pow(u[i][j],2) + pow(v[i][j],2));
    if(tmp <=  umax*1.0            ){a[i][j] = '9';}
    if(tmp <= (umax*0.9 + umin*0.1)){a[i][j] = '8';}
    if(tmp <= (umax*0.8 + umin*0.2)){a[i][j] = '7';}
    if(tmp <= (umax*0.7 + umin*0.3)){a[i][j] = '6';}
    if(tmp <= (umax*0.6 + umin*0.4)){a[i][j] = '5';}
    if(tmp <= (umax*0.5 + umin*0.5)){a[i][j] = '4';}
    if(tmp <= (umax*0.4 + umin*0.6)){a[i][j] = '3';}
    if(tmp <= (umax*0.3 + umin*0.7)){a[i][j] = '2';}
    if(tmp <= (umax*0.2 + umin*0.8)){a[i][j] = '1';}
    if(tmp <= (umax*0.1 + umin*0.9)){a[i][j] = '0';}
  } }

  printf("Re = %4.2f, ut = %5.3f, tau = %5.3f\n", re, ut, tau);

  for(j = ny - 1; j >= 1; j = j - 4){ for(i = 1; i <= nx - 1; i = i + 4){
    printf("%c",a[i][j]);
    }
    printf("\n");
  }

  //vorticity
 
  for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){
    omega[i][j] = 0.0;
  } }
  for(i = 1; i < nx; i++){ for(j = 1; j < ny; j++){
    omega[i][j] = (v[i+1][j  ] - v[i-1][j  ])*0.5
                - (u[i  ][j+1] - u[i  ][j-1])*0.5;
  } }

  sprintf(filename, "step%04d.vti", time);

  fp = fopen(filename, "w");
    fprintf(fp, "<?xml version=\"1.0\"?>\n");
    fprintf(fp, "<VTKFile type=\"ImageData\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
    fprintf(fp, "<ImageData WholeExtent=\"0 %d 0 %d 0 0\" Origin=\"0 0 0\" Spacing=\"1 1 1\">\n",nx-2, ny-2);
    fprintf(fp, "<Piece Extent=\"0 %d 0 %d 0 0\">\n", nx-2, ny-2);
    fprintf(fp, "<PointData Vectors=\"velocity\" Scalars=\"vorticity\">\n");

    fprintf(fp, "<DataArray type=\"Float64\" Name=\"velocity\" NumberOfComponents=\"3\" format=\"ascii\">\n");
    for (int j = 1; j < ny; j++) {
        for (int i = 1; i < nx; i++) {
            fprintf(fp, "%f %f 0.0\n", u[i][j], v[i][j]);
        }
    }
    fprintf(fp, "</DataArray>\n");

    fprintf(fp, "<DataArray type=\"Float64\" Name=\"vorticity\" NumberOfComponents=\"1\" format=\"ascii\">\n");
    for (int j = 1; j < ny; j++) {
        for (int i = 1; i < nx; i++) {
            fprintf(fp, "%f\n", omega[i][j]);
        }
    }
    fprintf(fp, "</DataArray>\n");

    fprintf(fp, "</PointData>\n");
    fprintf(fp, "<CellData>\n</CellData>\n");
    fprintf(fp, "</Piece>\n");
    fprintf(fp, "</ImageData>\n");
    fprintf(fp, "</VTKFile>\n");

    fclose(fp);

  fprintf(fp_pvd, "  <DataSet timestep=\"%d\" file=\"step%04d.vti\"/>\n", time, time);

  if(norm < 0.000000000001 && time > 10000){ exit(0); }

  } //loop1

  fprintf(fp_pvd, "  </Collection>\n</VTKFile>\n");
  fclose(fp_pvd);

  return 0;
}