//
#include "omp.h"
#include<iostream>
#include<cmath>
#include<cstdlib>
#include<iomanip>
#include<fstream>
#include<sstream>
#include<string>
#include<time.h>
//命名空间
using namespace std;
//常量定义
const int Q=9;
const int NX=64;
const int NY=64;
const double ERROR=1e-8;
const int Step_scr=1000;
const int Step_file=5000;
const int NXm=NX-1;
const int NYm=NY-1;
const double Cs=0.01; //Smagorinsky常数
const double Delta=1.0; //滤波宽度
const double Th=20.0;
const double Tl=1.0;
const double G=-9.8;
const double Ma=0.1;
const double Pr=0.71;
const double Ra=1e3;
#define BETA Ma*Ma*cscs/G/(Th-Tl)/Ly
//变量定义
int e[Q][2]={{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};
double w[Q]={4.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/36,1.0/36,1.0/36,1.0/36};
double u[NX+1][NY+1][2],u0[NX+1][NY+1][2];
double fs0[NX+1][NY+1][Q],fs[NX+1][NY+1][Q],fc[NX+1][NY+1][Q],feq[NX+1][NY+1][Q];
double gs0[NX+1][NY+1][Q],gs[NX+1][NY+1][Q],gc[NX+1][NY+1][Q],geq[NX+1][NY+1][Q];
double rho[NX+1][NY+1],t[NX+1][NY+1],t0[NX+1][NY+1],psi[NX+1][NY+1],ksi[NX+1][NY+1],p[NX+1][NY+1];
int i,j,k,ip,jp,n,ic,jc;
double c,cs,cscs,Re,dx,dy,Lx,Ly,dt,rho0,p0,tau_0,tau_1,t_0,niu,error1,error2,_u;
double psi1,psi2;
double tau_LES[NX+1][NY+1];
double nu[NY+1],nu_ave,nu_ave0,nu_sum; //热壁努塞尔数
time_t time1,time2;
//函数声明
void Init(); //done
void tau_calc();
void Stream(); //done
void Macro_calc(); //done
void Feq_calc(); //done
void Collide(); //done
void Edges(); //done
void Var_calc(); //计算需要额外输出至文件的变量
void Psi_calc(); //done
void Ksi_calc(); //done
void Nu_calc();
void Error(); //done
void Out_scr(); //done
void Out_file(int m); //done
int main()
{
omp_set_num_threads(2);
Init(); //初始化
time(&time2);
for (n=0;;n++)
{
Nu_calc();
ostringstream name5;
name5<<"Nu.dat";
ofstream out5(name5.str().c_str(),ios::app);
out5<<n<<" "<<nu_ave<<endl;
ostringstream mpoint;
mpoint<<"mpoint.dat";
ofstream mout(mpoint.str().c_str(),ios::app);
mout<<n<<" "<<u[int(2.0/16*NX)][int(13.0/16*NX)][0]/_u<<" "
<<u[int(2.0/16*NX)][int(2.0/16*NX)][1]/_u<<" "<<t[int(2.0/16*NX)][int(2.0/16*NX)]<<endl;
Stream(); //迁移
Macro_calc(); //宏观量统计
Feq_calc(); //平衡态分布函数
//tau_calc(); //大涡模拟
Collide(); //碰撞
Edges(); //边界条件
if (n%Step_scr==0)
{
Error();
nu_ave0=nu_ave;
Out_scr(); //屏幕显示输出
if (n>=Step_file && n%Step_file==0)
{
cout<<endl<<error1<<endl;
Var_calc(); //计算需要额外输出至文件的变量（流函数、涡量等）
Out_file(n); //数据文件输出
time1=time2;
time(&time2);
cout<<endl<<time2-time1<<" seconds cost in the past "<<Step_file<<" steps!"<<endl<<endl;
//cout<<"Change of the minimum value of the stream function: "
//<<setiosflags(ios::scientific)<<psi1-psi2<<endl<<"Center: "<<double(ic)/NX<<" , "<<double(jc)/NY<<endl<<endl;
if (error1<ERROR && error2<ERROR)
break;//收敛判断
//if (fabs(psi1-psi2)<1e-5)
//if (n>=1000000)
//break;
}
}
}
return 0;
}
void Init()//初始化
{
nu_sum=0;
nu_ave=1.0;
nu_ave0=0.0;
t_0=(Th+Tl)/2;
dx=1.0;
dy=1.0;
Lx=dx*double(NX);
Ly=dy*double(NY);
dt=dx;
c=dx/dt; //c设定为1
cscs=c*c/3.0;
cs=sqrt(cscs);
rho0=3.0;
tau_0=Ma/cs/dt*Lx*sqrt(Pr/Ra)+0.5;
tau_1=(tau_0-0.5)/Pr+0.5;
_u=cscs*dt*(tau_0-0.5)/Pr/Ly; //参考速度
for (i=0;i<=NX;i++)
for (j=0;j<=NY;j++)
{
u[i][j][0]=0;
u[i][j][1]=0;
u0[i][j][0]=0.1;
u0[i][j][1]=0.1;
t0[i][j]=0;
rho[i][j]=rho0;
t[i][j]=t_0;
t[0][j]=Th;
t[NX][j]=Tl;
psi[i][j]=0;
ksi[i][j]=0;
tau_LES[i][j]=0;
}
Feq_calc();
for (i=0;i<=NX;i++)
for (j=0;j<=NY;j++)
for (k=0;k<Q;k++)
{
fs0[i][j][k]=feq[i][j][k];
fs[i][j][k]=fs0[i][j][k];
fc[i][j][k]=fs[i][j][k];
gs0[i][j][k]=geq[i][j][k];
gs[i][j][k]=gs0[i][j][k];
gc[i][j][k]=gs[i][j][k];
}
cout<<"Initialiation Complete!"<<endl
<<"tau_0= "<<tau_0<<" tau_1= "<<tau_1<<endl;
ostringstream name5;
name5<<"Nu.dat";
ofstream out5(name5.str().c_str());
out5<<"VARIABLES=n,Nu\n";
{
ostringstream mpoint;
mpoint<<"mpoint.dat";
ofstream mout(mpoint.str().c_str());
mout<<"VARIABLES=n,U,V,t\n";
}
}
void Macro_calc()//计算宏观量
{
#pragma omp parallel for private(i,j,k)
for (i=1;i<NX;i++)
for (j=1;j<NY;j++)
{
rho[i][j]=0;
u[i][j][0]=0;
u[i][j][1]=0;
t[i][j]=0;
for(k=0;k<Q;k++)
{
rho[i][j]+=fs[i][j][k];
u[i][j][0]+=e[k][0]*fs[i][j][k];
u[i][j][1]+=e[k][1]*fs[i][j][k];
t[i][j]+=gs[i][j][k];
}
u[i][j][0]/=rho[i][j];
u[i][j][1]/=rho[i][j];
}
for (i=1;i<NX;i++)
{
rho[i][0]=rho[i][1];
rho[i][NY]=rho[i][NYm];
t[i][0]=t[i][1];
t[i][NY]=t[i][NYm];
}
for (j=0;j<=NY;j++)
{
rho[0][j]=rho[1][j];
rho[NX][j]=rho[NXm][j];
t[0][j]=Th;
t[NX][j]=Tl;
}
}
void Feq_calc()//计算平衡态分布函数
{
double eu,uu;
#pragma omp parallel for private(i,j,k,eu,uu)
for (i=0;i<=NX;i++)
for (j=0;j<=NY;j++)
{
uu=(u[i][j][0]*u[i][j][0]+u[i][j][1]*u[i][j][1]);
for (k=0;k<Q;k++)
{
eu=(e[k][0]*u[i][j][0]+e[k][1]*u[i][j][1]);
feq[i][j][k]=w[k]*rho[i][j]*(1.0+3.0*eu+4.5*eu*eu-1.5*uu);
geq[i][j][k]=w[k]*t[i][j]*(1.0+3.0*eu+4.5*eu*eu-1.5*uu);
}
}
}
void Stream()//迁移过程
{
#pragma omp parallel for private(i,j,k,ip,jp)
for (i=1;i<NX;i++)
for (j=1;j<NY;j++)
{
for (k=0;k<Q;k++)
{
ip=i-e[k][0];
jp=j-e[k][1];
fs[i][j][k]=fc[ip][jp][k];
gs[i][j][k]=gc[ip][jp][k];
}
}
#pragma omp parallel for private(i,j,k)
for (i=0;i<=NX;i++)
for (j=0;j<=NY;j++)
for (k=0;k<Q;k++)
{
fs0[i][j][k]=fc[i][j][k];
gs0[i][j][k]=gc[i][j][k];
}
}
void Var_calc()//计算需要额外输出至文件的变量
{
Psi_calc();
Ksi_calc();
Nu_calc();
for (i=0;i<=NX;i++)
for (j=0;j<=NY;j++)
p[i][j]=rho[i][j]*cscs;
}
void Psi_calc()//计算流函数
{
int t;
psi1=psi2;
psi2=0;
for (j=1;j<NY;j++)
{
psi[0][j]=0;
for (i=1;i<NX;i++)
{
psi[i][j]=0;
for (t=1;t<i;t++,t++)
psi[i][j]-=(u[t-1][j][1]+4.0*u[t][j][1]+u[t+1][j][1])/3.0*dx/_u/Lx;
if (t==i)
psi[i][j]-=(u[t-1][j][1]+u[t][j][1])/2.0*dx/_u/Lx;
if (psi[i][j]<psi2)
{
psi2=psi[i][j];
ic=i;jc=j;
}
}
}
}
void Ksi_calc()//计算涡量
{
for (i=1;i<NX;i++)
for (j=1;j<NY;j++)
ksi[i][j]=((u[i+1][j][1]-u[i-1][j][1])/dx-(u[i][j+1][0]-u[i][j-1][0])/dy)/2/_u*Lx;
}
void Collide()//碰撞
{
double Gy,eu,euy;
#pragma omp parallel for private(i,j,k)
for (i=1;i<NX;i++)
for (j=1;j<NY;j++)
for (k=0;k<Q;k++)
{
fc[i][j][k]=fs[i][j][k]-(fs[i][j][k]-feq[i][j][k])/tau_0;
Gy=rho[i][j]*BETA*(t[i][j]-t_0)*G;
//eu=e[k][0]*u[i][j][0]+e[k][1]*u[i][j][1];
//euy=(e[k][1]-u[i][j][1])/cscs+eu/cscs/cscs*e[k][1];
euy=(e[k][1]-u[i][j][1])/rho[i][j]/cscs;
//if (k==2||k==4)
fc[i][j][k]+=Gy*euy*fs[i][j][k];
gc[i][j][k]=gs[i][j][k]-(gs[i][j][k]-geq[i][j][k])/tau_1;
}
}
void Edges()//边界条件
{
for (j=1;j<NY;j++)//左右
for (k=0;k<Q;k++)
{
fc[0][j][k]=feq[0][j][k]+(fc[1][j][k]-feq[1][j][k]);
fc[NX][j][k]=feq[NX][j][k]+(fc[NXm][j][k]-feq[NXm][j][k]);
gc[0][j][k]=geq[0][j][k]+(gc[1][j][k]-geq[1][j][k]);
gc[NX][j][k]=geq[NX][j][k]+(gc[NXm][j][k]-geq[NXm][j][k]);
}
for (i=1;i<NX;i++)//上下
for (k=0;k<Q;k++)
{
fc[i][0][k]=feq[i][0][k]+(fc[i][1][k]-feq[i][1][k]);
fc[i][NY][k]=feq[i][NY][k]+(fc[i][NYm][k]-feq[i][NYm][k]);
gc[i][0][k]=geq[i][0][k]+(gc[i][1][k]-geq[i][1][k]);
gc[i][NY][k]=geq[i][NY][k]+(gc[i][NYm][k]-geq[i][NYm][k]);
}
for (k=0;k<Q;k++)//角点
{
fc[0][0][k]=feq[0][0][k]+(fc[1][1][k]-feq[1][1][k]);
fc[0][NY][k]=feq[0][NY][k]+(fc[1][NYm][k]-feq[1][NYm][k]);
fc[NX][0][k]=feq
