// cavity flow (Cumulant)
// 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
// f  : distribution function
// rm : raw moment
// cm : central moment
// cu : cumurant
// nx : number of grid points (x-axis)
// ny : number of grid points (y-axis)
// 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], ftmp[9][DIM][DIM], cx[9], cy[9];
  double rm[9][DIM][DIM], cm[9][DIM][DIM], cu[9][DIM][DIM];
  double ut = 0.050, vt = 0.0, umax, umin, tmp, u2, mu, norm, tau, re = 10000;
  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++){ for(k = 0; k <= 8; k++){
    f[k][i][j] = 0.0;
  } } }

// calculation start
  for(loop1 = 0; loop1 < 30; loop1++){
  for(loop2 = 0; loop2 < 10000; 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];
    } }

// cumulant (distribution function - > raw moment)
//  rm[0]:m_00, rm[1]:m_10, rm[2]:m_01
//  rm[3]:m_11, rm[4]:m_20, rm[5]:m_02
//  rm[6]:m_21, rm[7]:m_12, rm[8]:m_22

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

      for( k = 0; k <= 8; k++){
        rm[0][i][j] += f[k][i][j];
        rm[1][i][j] += f[k][i][j]*cx[k];
        rm[2][i][j] += f[k][i][j]*cy[k];
        rm[3][i][j] += f[k][i][j]*cx[k]      *cy[k];
        rm[4][i][j] += f[k][i][j]*cx[k]*cx[k];
        rm[5][i][j] += f[k][i][j]            *cy[k]*cy[k];
        rm[6][i][j] += f[k][i][j]*cx[k]*cx[k]*cy[k];
        rm[7][i][j] += f[k][i][j]*cx[k]      *cy[k]*cy[k];
        rm[8][i][j] += f[k][i][j]*cx[k]*cx[k]*cy[k]*cy[k];
      } 
        rm[0][i][j] += 1.0;
        rm[4][i][j] += 1.0/3.0;
        rm[5][i][j] += 1.0/3.0;
        rm[8][i][j] += 1.0/9.0;
    } }

// cumulant (raw moment - > central moment)
//  cm[0]:kappa_00, cm[1]:kappa_10, cm[2]:kappa_01
//  cm[3]:kappa_11, cm[4]:kappa_20, cm[5]:kappa_02
//  cm[6]:kappa_21, cm[7]:kappa_12, cm[8]:kappa_22
    for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){
      for( k = 0; k <= 8; k++){
        cm[k][i][j] = 0.0;
      }
        cm[0][i][j] = rm[0][i][j];
        cm[1][i][j] = 0.0;
        cm[2][i][j] = 0.0;
        cm[3][i][j] =-rm[0][i][j]*u[i][j]*v[i][j] + rm[3][i][j];
        cm[4][i][j] =-rm[0][i][j]*u[i][j]*u[i][j] + rm[4][i][j];
        cm[5][i][j] =-rm[0][i][j]*v[i][j]*v[i][j] + rm[5][i][j];
        cm[6][i][j] = 2.0*rm[0][i][j]*u[i][j]*u[i][j]*v[i][j] - 2.0*rm[3][i][j]*u[i][j] - rm[4][i][j]*v[i][j] + rm[6][i][j];
        cm[7][i][j] = 2.0*rm[0][i][j]*u[i][j]*v[i][j]*v[i][j] - 2.0*rm[3][i][j]*v[i][j] - rm[5][i][j]*u[i][j] + rm[7][i][j];
        cm[8][i][j] =-3.0*rm[0][i][j]*u[i][j]*u[i][j]*v[i][j]*v[i][j] + rm[5][i][j]*u[i][j]*u[i][j] + rm[4][i][j]*v[i][j]*v[i][j]
                    + 4.0*rm[3][i][j]*u[i][j]*v[i][j] - 2.0*rm[7][i][j]*u[i][j] - 2.0*rm[6][i][j]*v[i][j] + rm[8][i][j];
    } }

// cumulant (central moment - > cumulant + collision)
//  cu[0]:cumulant_00, cu[1]:cumulant_10, cu[2]:cumulant_01
//  cu[3]:cumulant_11, cu[4]:cumulant_20, cu[5]:cumulant_02
//  cu[6]:cumulant_21, cu[7]:cumulant_12, cu[8]:cumulant_22
    for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){
      for( k = 0; k <= 8; k++){
        cu[k][i][j] = 0.0;
      }
        cu[3][i][j] =  cm[3][i][j] - cm[3][i][j]/tau;
        cu[4][i][j] = (cm[4][i][j] - cm[5][i][j]) - (cm[4][i][j] - cm[5][i][j])/tau;
        cu[5][i][j] = rm[0][i][j]*2.0/3.0;
        cu[6][i][j] = 0.0;
        cu[7][i][j] = 0.0;
        cu[8][i][j] = 0.0;
    } }

// cumulant (cumulant - > central moment)
//  cm[0]:kappa_00, cm[1]:kappa_10, cm[2]:kappa_01
//  cm[3]:kappa_11, cm[4]:kappa_20, cm[5]:kappa_02
//  cm[6]:kappa_21, cm[7]:kappa_12, cm[8]:kappa_22
    for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){
        cm[3][i][j] =  cu[3][i][j];
        cm[4][i][j] =( cu[4][i][j] + cu[5][i][j])*0.5;
        cm[5][i][j] =(-cu[4][i][j] + cu[5][i][j])*0.5;
        cm[6][i][j] =  cu[6][i][j];
        cm[7][i][j] =  cu[7][i][j];
        cm[8][i][j] =  cu[8][i][j] + cm[4][i][j]*cm[5][i][j]/rm[0][i][j] + 2.0 *cm[3][i][j]*cm[3][i][j]/rm[0][i][j];
    } }

// cumulant (central moment -> raw moment)
//  rm[0]:m_00, rm[1]:m_10, rm[2]:m_01
//  rm[3]:m_11, rm[4]:m_20, rm[5]:m_02
//  rm[6]:m_21, rm[7]:m_12, rm[8]:m_22
    for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){
        rm[1][i][j] = rm[0][i][j]*u[i][j];
        rm[2][i][j] = rm[0][i][j]*v[i][j];
        rm[3][i][j] = cm[3][i][j] + rm[0][i][j]*u[i][j]*v[i][j];
        rm[4][i][j] = cm[4][i][j] + rm[0][i][j]*u[i][j]*u[i][j];
        rm[5][i][j] = cm[5][i][j] + rm[0][i][j]*v[i][j]*v[i][j];
        rm[6][i][j] = 2.0*cm[3][i][j]*u[i][j] + rm[0][i][j]*u[i][j]*u[i][j]*v[i][j] + cm[6][i][j] + cm[4][i][j]*v[i][j];
        rm[7][i][j] = 2.0*cm[3][i][j]*v[i][j] + rm[0][i][j]*u[i][j]*v[i][j]*v[i][j] + cm[7][i][j] + cm[5][i][j]*u[i][j];
        rm[8][i][j] = rm[0][i][j]*u[i][j]*u[i][j]*v[i][j]*v[i][j] + cm[5][i][j]*u[i][j]*u[i][j] + cm[4][i][j]*v[i][j]*v[i][j]
        + 4.0*cm[3][i][j]*u[i][j]*v[i][j] + 2.0*cm[7][i][j]*u[i][j] + 2.0*cm[6][i][j]*v[i][j] + cm[8][i][j];
    } }

// cumulant  (raw moment  - >  distribution function)
//  rm[0]:m_00, rm[1]:m_10, rm[2]:m_01
//  rm[3]:m_11, rm[4]:m_20, rm[5]:m_02
//  rm[6]:m_21, rm[7]:m_12, rm[8]:m_22
    for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){
        f[0][i][j] = (rm[0][i][j] - 1.0) - (rm[5][i][j] - 1.0/3.0) - (rm[4][i][j] - 1.0/3.0) + (rm[8][i][j] - 1.0/9.0);
        f[1][i][j] = ((rm[5][i][j] - 1.0/3.0) - (rm[8][i][j] - 1.0/9.0))*0.5 + ( rm[2][i][j] - rm[6][i][j])*0.5;
        f[2][i][j] = ((rm[5][i][j] - 1.0/3.0) - (rm[8][i][j] - 1.0/9.0))*0.5 - ( rm[2][i][j] - rm[6][i][j])*0.5;
        f[3][i][j] = ((rm[4][i][j] - 1.0/3.0) - (rm[8][i][j] - 1.0/9.0))*0.5 + (-rm[1][i][j] + rm[7][i][j])*0.5;
        f[4][i][j] = ((rm[4][i][j] - 1.0/3.0) - (rm[8][i][j] - 1.0/9.0))*0.5 - (-rm[1][i][j] + rm[7][i][j])*0.5;
        f[5][i][j] = (-rm[3][i][j] + rm[8][i][j] - 1.0/9.0)*0.25 + (-rm[7][i][j] + rm[6][i][j])*0.25;
        f[6][i][j] = ( rm[3][i][j] + rm[8][i][j] - 1.0/9.0)*0.25 + ( rm[7][i][j] + rm[6][i][j])*0.25;
        f[7][i][j] = ( rm[3][i][j] + rm[8][i][j] - 1.0/9.0)*0.25 - ( rm[7][i][j] + rm[6][i][j])*0.25;
        f[8][i][j] = (-rm[3][i][j] + rm[8][i][j] - 1.0/9.0)*0.25 - (-rm[7][i][j] + rm[6][i][j])*0.25;
    } }

    // 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];

      f[2][i][ny-1] = f[1][i ][ny] - (cx[1]*ut + cy[1]*vt)*2.0/3.0;
      f[7][i][ny-1] = f[6][ip][ny] - (cx[6]*ut + cy[6]*vt)/6.0;
      f[8][i][ny-1] = f[5][im][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] = (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] = (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] = (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] = (3.0*tmp +9.0/2.0*tmp*tmp -3.0/2.0*u2)/36.0;

    f[6][nx-1][   1] = 0.0; f[7][nx-1][   1] = 0.0;
    f[5][   1][   1] = 0.0; f[8][   1][   1] = 0.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] += f[k][i][j];
          u[i][j] += f[k][i][j]*cx[k];
          v[i][j] += f[k][i][j]*cy[k];
      } 
    } }

    for(i = 0; i <= nx; i++){ for(j = 0; j <= ny; j++){
        rho[i][j] += 1.0;
      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
  
  printf("Time = %d, Norm = %8.6e\n", time, norm);

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