The code below is posted as a support for the following thread on the forum.
/* FVLBM code Problem: Poiseuille flow Assumptions: 1- this model is cell vertex finite volume metohd 2- each vertex has 8 neighbors 3- dx is the same between all vertices 4 3 2 \ | / \|/ 5--0--1 /|\ / | \ 6 7 8 this is the numbering of neighbors and also the distributio functions Pnt is a help sruct to hold a pier of data */ #include <iostream> #include <fstream> #include <math.h> using namespace std; template <typename T>struct Pnt { Pnt(){v[0]=0,v[1]=0;} Pnt(T v1, T v2) { v[0]=v1,v[1]=v2;} T v[2]; T& operator[](int index){ return v[index];} }; void Q9(int& q,int**& c, double*& t, int*& ops, double& csq ) { q=9; t=new double[q]; ops=new int[q]; c=new int*[2]; for (int iD=0; iD<2; ++iD) c[iD]=new int[q]; t[0]=4./9.; t[1]=1./9.; t[2]=1./9.; t[3]=1./36.; t[4]=1./36; t[5]=1./9.; t[6]=1./9.; t[7]=1./36.; t[8]=1./36; c[0][0]=0; c[0][1]=-1; c[0][2]= 0; c[0][3]=-1; c[0][4]=-1; c[1][0]=0; c[1][1]= 0; c[1][2]=-1; c[1][3]= 1; c[1][4]=-1; c[0][5]=+1; c[0][6]= 0; c[0][7]=+1; c[0][8]=+1; c[1][5]= 0; c[1][6]=+1; c[1][7]=-1; c[1][8]=+1; ops[0]=0; ops[1]=5; ops[2]=6; ops[3]=7; ops[4]=8; ops[5]=1; ops[6]=2; ops[7]=3; ops[8]=4; csq=1./3.; } void Q9v(int& q,int**& c, double*& t, int*& ops, double& csq ) { q=9; t=new double[q]; ops=new int[q]; c=new int*[2]; for (int iD=0; iD<2; ++iD) c[iD]=new int[q]; t[0]=4./9.; t[1]=1./9.; t[2]=1./36.; t[3]=1./9.; t[4]=1./36; t[5]=1./9.; t[6]=1./36.; t[7]=1./9.; t[8]=1./36; c[0][0]=0; c[0][1]=+1; c[0][2]=+1; c[0][3]= 0; c[0][4]=-1; c[1][0]=0; c[1][1]= 0; c[1][2]=+1; c[1][3]= 1; c[1][4]=+1; c[0][5]=-1; c[0][6]=-1; c[0][7]= 0; c[0][8]=+1; c[1][5]= 0; c[1][6]=-1; c[1][7]=-1; c[1][8]=-1; ops[0]=0; ops[1]=5; ops[2]=6; ops[3]=7; ops[4]=8; ops[5]=1; ops[6]=2; ops[7]=3; ops[8]=4; csq=1./3.; } class Core_2D { private: int half ; int nextX ,nextY; double*** f; double*** ftemp; //int nX,nY; int q; int d; int** c; double* t; int* ops; double csq; double temp; ofstream out; double** Rho, **U [2]; double ( Core_2D::*Feq ) ( int iF, double rho, Pnt<double> u, double uSqr ); Pnt<double> ** P; public: double Omega; double Visc; int nX,nY,nT; Core_2D(int I, int J, int** ck,double* wk, int* opsit ,double& cs, int Q , int type) { out.open("out.txt",ios::trunc); d=2; q=Q; half = (Q-1)/2; nX=I; nY=J; csq=cs; f =new double** [q]; ftemp =new double** [q]; ops=new int [q]; for(int iF=0;iF<q;iF++) { f[iF]=new double* [nX]; ftemp[iF]=new double* [nX]; for(int iX=0;iX<nX;iX++) { f[iF][iX]=new double [nY]; ftemp[iF][iX]=new double [nY]; } } c=new int* [d]; for(int iD=0;iD<d;iD++) c[iD]=ck[iD]; t=new double [q]; for(int iF=0;iF<q;iF++) { t[iF]=wk[iF]; ops[iF]=opsit[iF]; } if ( type==1 ) Feq= &Core_2D::fEq1; else if ( type==2 ) Feq= &Core_2D::fEq2; else if ( type==4 ) Feq= &Core_2D::fEq4; Rho=new double* [nX]; U[0]=new double* [nX]; U[1]=new double* [nX]; P=new Pnt<double>* [nX]; for(int iX=0;iX<nX;iX++) { Rho[iX]=new double [nY]; U[0][iX]=new double [nY]; U[1][iX]=new double [nY]; P[iX]=new Pnt<double> [nY]; } } void Init_Geo(double L, double dx) { for(int iX=0;iX<nX;iX++) { for(int iY=0;iY<nY; iY++) { P[iX][iY][0]=iX*dx; P[iX][iY][1]=iY*dx; } } } // initual all the field with rho and u void Init (double rho, Pnt<double> u) { for (int iX=0; iX<nX; ++iX) for (int iY=0; iY<nY; ++iY) { for (int iF=0; iF<q; ++iF) f[iF][iX][iY]= (this->*Feq)(iF, rho, u,u[0]*u[0]+u[1]*u[1]); Rho[iX][iY]=rho; U[0][iX][iY]=u[0]; U[1][iX][iY]=u[1]; } } // initual all the field with rho void Init (double rho) { for (int iX=0; iX<nX; ++iX) for (int iY=0; iY<nY; ++iY) { for (int iF=0; iF<q; ++iF) f[iF][iX][iY]= rho*t[iF]; Rho[iX][iY]=rho; } } void Init (Pnt<int> ix, Pnt<int> iy, double rho, Pnt<double> u) { for (int iX=ix[0]; iX<=ix[1]; ++iX) for (int iY=iy[0]; iY<=iy[1]; ++iY) { for (int iF=0; iF<q; ++iF) f[iF][iX][iY]= (this->*Feq)(iF, rho, u,u[0]*u[0]+u[1]*u[1]); Rho[iX][iY]=rho; U[0][iX][iY]=u[0]; U[1][iX][iY]=u[1]; } } void Init (Pnt<int> ix, Pnt<int> iy, double rho) { for (int iX=ix[0]; iX<=ix[1]; ++iX) for (int iY=iy[0]; iY<=iy[1]; ++iY) Init (rho); } void computeRho(int iX, int iY, double& rho) { rho = 0.; for (int iF=0; iF<q; ++iF) { rho =rho +f[iF][iX][iY]; } } void computeU(int iX, int iY, double& rho, Pnt<double>& u, double& uSqr) { uSqr = 0.; for(int iD=0; iD<d; ++iD) { u[iD] = 0.; for (int iF=0; iF<q; ++iF) { u[iD] += f[iF][iX][iY] * c[iD][iF]; } //if(rho==0) // u[iD]=0; //else u[iD] /= rho; uSqr += u[iD]*u[iD]; } } double fEq1 (int iF, double rho, Pnt<double> u, double uSqr) { return rho*t[iF]; } double fEq2(int iF, double rho, Pnt<double> u, double uSqr) { double c_u = 0.; // scalar product between c_{iF} and u for (int iD=0; iD<d; iD++) { c_u += c[iD][iF] * u[iD]; } return rho*t[iF]*(1+c_u/csq); } double fEq4(int iF, double rho, Pnt<double> u, double uSqr) { double c_u = 0.; // scalar product between c_{iF} and u for (int iD=0; iD<d; iD++) { c_u += c[iD][iF] * u[iD]; } return rho*t[iF]*(1. + c_u/csq + c_u*c_u/csq/csq/2.0 - uSqr/csq/2.0); } void FVon(int Ix , int Iy, double dx, double dt, double tau, int nb_min, int nb_max) { const int node=8; bool isBoundary = (nb_max-nb_min==node)?false:true ; double Edge_Mid[node+1][2]; double SubVol_Cntr[node+1][2]; double f_Edge_Mid[node+1]; double f_SubVol_Cntr[node]; double feq_Edge_Mid[node+1]; double feq_SubVol_Cntr[node+1]; double Ynorm[9][2]; double Xnorm[9][2]; int N; int Nb; int Nb_min=(Nb_min+node-1)%node+1; int Nb_max=(Nb_max+node-1)%node+1;; Pnt<double> u; double uSqr; Edge_Mid [0 ][0]=P[Ix] [Iy][0]; Edge_Mid [0 ][1]=P[Ix] [Iy][1]; SubVol_Cntr[0 ][0]=P[Ix] [Iy][0]; SubVol_Cntr[0 ][1]=P[Ix] [Iy][1]; for(int iF=0; iF<q ;iF++) { f_Edge_Mid [0 ]=f[iF][Ix][Iy]; f_SubVol_Cntr[0 ]=f[iF][Ix][Iy]; u[0]=U[0][Ix][Iy];//+0.000012/omega; u[1]=U[1][Ix][Iy]; uSqr=u[0]*u[0]+u[1]*u[1]; feq_Edge_Mid [0 ]=(this->*Feq)(iF, Rho[Ix][Iy] ,u , uSqr); feq_SubVol_Cntr[0 ]=(this->*Feq)(iF, Rho[Ix][Iy] ,u , uSqr); Ynorm[0][0]= Edge_Mid [0 ][0] - Edge_Mid [nb_min][0]; Ynorm[0][1]= Edge_Mid [nb_max][0] - Edge_Mid [0 ][0]; Xnorm[0][0]= Edge_Mid [nb_min][1] - Edge_Mid [0 ][1]; Xnorm[0][1]= Edge_Mid [0 ][1] - Edge_Mid [nb_max][1]; for(int nb=nb_min; nb<nb_max ;nb++) { N =(nb+node-1)%node+1; Nb=(nb+node)%node+1; Edge_Mid [N ][0]=( P[Ix ] [Iy ][0]+ P[Ix+c[0][N ] ] [Iy+c[1][N ] ][0] ) /2.0; SubVol_Cntr[N ][0]=( P[Ix ] [Iy ][0]+ P[Ix+c[0][N ] ] [Iy+c[1][N ] ][0]+ P[Ix+c[0][Nb] ] [Iy+c[1][Nb] ][0] ) /3.0; Edge_Mid [N ][1]=( P[Ix ] [Iy ][1]+ P[Ix+c[0][N ] ] [Iy+c[1][N ] ][1] ) /2.0; SubVol_Cntr[N ][1]=( P[Ix ] [Iy ][1]+ P[Ix+c[0][N ] ] [Iy+c[1][N ] ][1]+ P[Ix+c[0][Nb] ] [Iy+c[1][Nb] ][1] ) /3.0; } for(int nb=nb_min; nb<nb_max ;nb++) { N =(nb+node-1)%node+1; Nb=(nb+node)%node+1; f_Edge_Mid[N ] =(f[iF][Ix][Iy]+f[iF][Ix+c[0][N ]][Iy+c[1][N ]])/2.0; f_SubVol_Cntr[N ]=( f[iF][Ix ][Iy ] + f[iF][Ix+c[0][N ] ][Iy+c[1][N ]] + f[iF][Ix+c[0][Nb] ][Iy+c[1][Nb]] )/3.0; u[0]=U[0][Ix][Iy];//+0.000012/omega; u[1]=U[1][Ix][Iy]; uSqr=u[0]*u[0]+u[1]*u[1]; double f0=(this->*Feq)(iF, Rho[Ix][Iy] ,u , uSqr); u[0]=U[0][Ix+c[0][N ]][Iy+c[1][N ]];//+0.000012/omega; u[1]=U[1][Ix+c[0][N ]][Iy+c[1][N ]]; uSqr=u[0]*u[0]+u[1]*u[1]; double f1=(this->*Feq)(iF, Rho[Ix+c[0][N ]][Iy+c[1][N ]] ,u , uSqr); u[0]=U[0][Ix+c[0][Nb]][Iy+c[1][Nb]];//+0.000012/omega; u[1]=U[1][Ix+c[0][Nb]][Iy+c[1][Nb]]; uSqr=u[0]*u[0]+u[1]*u[1]; double f2=(this->*Feq)(iF, Rho[Ix+c[0][Nb]][Iy+c[1][Nb]] ,u , uSqr); feq_Edge_Mid[N ] =(f0+f1)/2.0; feq_SubVol_Cntr[N ]=( f0 + f1 + f2 )/3.0; } //double Ynorm=xi - xi+1 //double Xnorm=yi+1 - yi for(int nb=nb_min; nb<nb_max ;nb++) { N =(nb+node-1)%node+1; Nb=(nb+node)%node+1; Ynorm[N ][0]= Edge_Mid [N ][0] - SubVol_Cntr[N ][0]; Ynorm[N ][1]=-Edge_Mid [Nb][0] + SubVol_Cntr[N ][0]; Xnorm[N ][0]= SubVol_Cntr[N ][1] - Edge_Mid [N ][1] ; Xnorm[N ][1]= Edge_Mid [Nb][1] - SubVol_Cntr[N ][1]; } double flux=0; for(int nb=nb_min; nb<nb_max ;nb++) { N =(nb+node-1)%node+1; Nb=(nb+node)%node+1; flux+=c[0][iF] * Xnorm[N ][0] * (f_Edge_Mid[N ]+ f_SubVol_Cntr[N ])/2.0; //x flux+=c[1][iF] * Ynorm[N ][0] * (f_Edge_Mid[N ]+ f_SubVol_Cntr[N ])/2.0; //y flux+=c[0][iF] * Xnorm[N ][1] * (f_Edge_Mid[Nb]+ f_SubVol_Cntr[N ])/2.0; //x flux+=c[1][iF] * Ynorm[N ][1] * (f_Edge_Mid[Nb]+ f_SubVol_Cntr[N ])/2.0; //y } if(isBoundary) { Nb=(nb_min+node)%node+1; flux+=c[0][iF] * Xnorm[0][0] * (f_Edge_Mid[Nb]+ f_SubVol_Cntr[0])/2.0; //x flux+=c[1][iF] * Ynorm[0][0] * (f_Edge_Mid[Nb]+ f_SubVol_Cntr[0])/2.0; //y Nb=(nb_max+node-1)%node+1; flux+=c[0][iF] * Xnorm[0][1] * (f_Edge_Mid[Nb]+ f_SubVol_Cntr[0])/2.0; //x flux+=c[1][iF] * Ynorm[0][1] * (f_Edge_Mid[Nb]+ f_SubVol_Cntr[0])/2.0; //y } double collision=0; for(int nb=nb_min; nb<nb_max ;nb++) { N =(nb+node-1)%node+1; Nb=(nb+node)%node+1; collision+= -dx/12.0/tau*(f_Edge_Mid [0 ]+ f_Edge_Mid [N ]+ f_SubVol_Cntr [N ] -feq_Edge_Mid [0 ]- feq_Edge_Mid[N ]- feq_SubVol_Cntr[0 ])/3.0; collision+= -dx*dx/12.0/tau*(f_Edge_Mid [0 ]+ f_Edge_Mid [Nb]+ f_SubVol_Cntr [N ] -feq_Edge_Mid [0 ]- feq_Edge_Mid[Nb]- feq_SubVol_Cntr[0 ])/3.0; } ftemp[iF][Ix][Iy]=f[iF][Ix][Iy]+dt/dx/dx*(collision-flux)+c[0][iF]/6.0*0.000012; } } void FVon(Pnt<int> ix, Pnt<int> iy, double dx, double dt, double tau) { for (int iX=ix[0]; iX<=ix[1]; ++iX) for (int iY=iy[0]; iY<=iy[1]; ++iY) FVon(iX , iY, dx, dt, tau,1,9); //for (int iX=ix[0]; iX<=ix[1]; ++iX) for (int iY=iy[0]; iY<=iy[1]; ++iY) //w FVon(0 , iY, dx, dt, tau,7,8+3); //for (int iX=ix[0]; iX<=ix[1]; ++iX) for (int iY=iy[0]; iY<=iy[1]; ++iY) //e FVon(nX-1 , iY, dx, dt, tau,3,7); for (int iX=ix[0]; iX<=ix[1]; ++iX) //s //for (int iY=iy[0]; iY<=iy[1]; ++iY) FVon(iX , 0, dx, dt, tau,1,5); for (int iX=ix[0]; iX<=ix[1]; ++iX) //n //for (int iY=iy[0]; iY<=iy[1]; ++iY) FVon(iX , nY-1, dx, dt, tau,5,8); for (int iX=ix[0]; iX<=ix[1]; ++iX) for (int iY=iy[0]; iY<=iy[1]; ++iY) for(int iF=0; iF<9 ;iF++) f[iF][iX][iY]=ftemp[iF][iX][iY]; } void collosion(double omega) { double uSqr; Pnt<double> u; for (int iX=0; iX<nX; ++iX) { for (int iY=0; iY<nY; ++iY) { u[0]=U[0][iX][iY]+0.000012/omega; //2.604e-5/omega; u[1]=U[1][iX][iY]; uSqr=U[0][iX][iY]*U[0][iX][iY]+U[1][iX][iY]*U[1][iX][iY]; uSqr=u[0]*u[0]+u[1]*u[1]; for (int iF=0; iF<q; ++iF) { f[iF][iX][iY] *= (1.-omega); f[iF][iX][iY] += omega * (this->*Feq)(iF, Rho[iX][iY], u, uSqr); } for (int iF=1; iF<q; ++iF) { if( c[0][iF]==-1 || (c[0][iF]==0 && c[1][iF] !=1) ) swap(f[iF][iX][iY], f[ops[iF]][iX][iY]); } } } } void streaming() { for (int iX=0; iX<nX; ++iX) { for (int iY=0; iY<nY; ++iY) { for (int iF=1; iF<q; ++iF) { nextX = iX + c[0][iF]; nextY = iY + c[1][iF]; if( c[0][iF]==-1 || (c[0][iF]==0 && c[1][iF] !=1) ) if (nextX>=0 && nextX<nX && nextY>=0 && nextY<nY) { swap(f[ops[iF]][iX][iY], f[iF][nextX][nextY]); } } } } } void update(int iX, int iY) { double ssum=0,usum=0,vsum=0; for (int iF=0; iF<q; ++iF) { ssum+=f[iF][iX][iY]; usum+=f[iF][iX][iY]*c[0][iF]; vsum+=f[iF][iX][iY]*c[1][iF]; } Rho[iX][iY]=ssum; if (ssum==0) { U[0][iX][iY]=0; U[1][iX][iY]=0; } else { U[0][iX][iY]=usum/ssum; U[1][iX][iY]=vsum/ssum; } } void update(Pnt<int> ix, Pnt<int> iy) { for (int iX=ix[0]; iX<=ix[1]; ++iX) for (int iY=iy[0]; iY<=iy[1]; ++iY) update(iX,iY); } void update() { for (int iX=0; iX<nX; ++iX) for (int iY=0; iY<nY; ++iY) update(iX,iY); } void BounceBackBoundOn(int iX , int iY, Pnt<int> IF) { for (int iF=1; iF<q; ++iF) //if( ( c[0][iF]==IF[0] && c[0][iF]!=0 ) || ( c[1][iF]==IF[1] && c[1][iF]!=0 ) ) if( Dir(iF, IF) ) f[iF][iX][iY]=f[ops[iF]][iX][iY]; update(iX,iY); } void BounceBackBoundOn(Pnt<int> ix, Pnt<int> iy, Pnt<int> IF) { for (int iX=ix[0]; iX<=ix[1]; ++iX) for (int iY=iy[0]; iY<=iy[1]; ++iY) BounceBackBoundOn(iX ,iY, IF); } void OpenBoundOn(int iX , int iY, Pnt<int> IF) { //for (int iF=0 ;iF<q; ++iF) // if( (c[0][iF]==-IF[0] && c[0][iF]!=0 ) || ( c[1][iF]==-IF[1] && c[1][iF]!=0 ) ) //f[iF][iX][iY]=4.0/3.0*f [iF] [ iX+IF[0] ] [ iY+IF[1] ] - 1.0/3.0*f[iF] [ iX+2*IF[0] ] [ iY+2*IF[1]]; //f[iF][iX][iY]=f[iF][iX+IF[0]][iY+IF[1]]; for (int iF=0 ;iF<q; ++iF) if( (c[0][iF]==IF[0] && c[0][iF]!=0 ) || ( c[1][iF]==IF[1] && c[1][iF]!=0 ) ) //f[iF][iX][iY]=4.0/3.0*f [iF] [ iX+IF[0] ] [ iY+IF[1] ] - 1.0/3.0*f[iF] [ iX+2*IF[0] ] [ iY+2*IF[1]]; f[iF][iX][iY]=f[iF][iX+IF[0]][iY+IF[1]]; update(iX,iY); for (int iD=0 ;iD<d; ++iD) if(IF[iD]==0) U[iD][iX][iY]=0; } void OpenBoundOn(Pnt<int> ix, Pnt<int> iy, Pnt<int> IF) { for (int iX=ix[0]; iX<=ix[1]; ++iX) for (int iY=iy[0]; iY<=iy[1]; ++iY) OpenBoundOn(iX ,iY, IF); } void TecPrint(int t) { double rho; Pnt<double> u; //ofstream out ( "out.txt",ios::ate ); double** str; str=new double*[nX]; for(int i=0; i<nX ; i++) str[i]=new double [nY]; for(int i=0; i<nX ; i++) for(int j=0; j<nY ; j++) str[i][j]=0; for(int i=0; i<nX ; i++) for(int j=0; j<nY ; j++) str[i][j]=IntgU( 0 , i , j) - IntgV(0 , 0 , i); out<<" VARIABLES = \"x\" \"y\" \"r\" \"u\" \"v\" \"s\" "<<endl; out<<" ZONE T= \"Z"<< t<<"\" "<<endl; out<<"I="<<nY<<", J="<<nX<<", K="<<1<<endl; for ( int iX=0;iX<nX;iX++ ) { for ( int iY=0;iY<nY;iY++ ) { computeRho(iX,iY,rho); computeU(iX, iY,rho, u, temp); out<<iX<<"\t"<<iY<<"\t"<<rho <<"\t"<<u[0]<<"\t"<<u[1]<<"\t"<<str[iX][iY] <<endl; } out<<endl; } //system ( "wgnuplot 3plot.plt -" ); } bool Dir(int iF, Pnt<int> IF) { if(IF[0]*IF[1]) return ( c[0][iF]!=-IF[0] ) && ( c[1][iF]!=-IF[1]) ; else return (c[0][iF]==IF[0] && c[0][iF]!=0 ) || ( c[1][iF]==IF[1] && c[1][iF]!= 0 ); //return (c2[0][iF]==IF[0] && c2[0][iF]!=-IF[0] ) || ( c2[1][iF]==IF[1] && c2[1][iF]!=-IF[1]); } double IntgV(int X0 , int Y0 , int X) { double I=0; for(int i=X0 ; i<X ; i++) I+= (U[1][i][Y0]+U[1][i+1][Y0])/2.0; return I; } double IntgU(int Y0 , int X , int Y) { double I=0; for(int j=Y0 ; j<Y ; j++) I+= (U[0][X][j]+U[0][X][j+1])/2.0; return I; } }; int main() { int q2; double* t2; int* ops2; int** c2; double csq2; int nX=32; int nY=nX; Q9v(q2,c2,t2,ops2,csq2 ); double dx_lu=0.5; double dt_lu=0.001; double tau_lu=0.1; double visc=tau_lu/3.0; //double omega2=1./(visc/csq2+0.5); //for FDLBM double omega2=1/tau_lu; //for FVLBM cout<< 1/omega2<<" "<<endl; char syspause; cin>>syspause; Core_2D f2(nX,nY ,c2,t2,ops2,csq2,q2,4); f2.Init_Geo(nX,dx_lu); //init the field with rho=1 and u=0 f2.Init (1, Pnt<double>(0.0,0.00) ); for(int i=0;i<6000;i++) { cout<<i<<endl; //f2.collosion(omega2); //for FDLBM //f2.streaming(); //for FDLBM f2.FVon(Pnt<int>(1,nX-2), Pnt<int>(1,nY-2), dx_lu, dt_lu, tau_lu); //for FVLBM f2.update(); f2.OpenBoundOn (Pnt<int>(0,0) , Pnt<int>(0,nY-1) , Pnt<int>(1,0)); //w f2.BounceBackBoundOn(Pnt<int>(0,nX-1) , Pnt<int>(nY-1,nY-1) , Pnt<int>(0,-1)); //n f2.BounceBackBoundOn(Pnt<int>(0,nX-1) , Pnt<int>(0,0) , Pnt<int>(0,+1)); //s f2.OpenBoundOn (Pnt<int>(nX-1,nX-1) , Pnt<int>(0,nY-1) , Pnt<int>(-1,0)); //e if(i%500==0) f2.TecPrint(i); } f2.TecPrint(400); }