// 方腔自然对流.cpp : Defines the entry point for the console application.
//不可压热格子模型——基于Boussinesq假设的耦合双分布函数模型


#include "stdafx.h"
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>

using namespace std;

const int Q = 9;//选择D2Q9模型
const int NX = 100;//x方向的格子数
const int NY = 100;//y方向的格子数
const int TH = 1;//左边界的壁面温度
const int TL = 0;//右边界的壁面温度
const double Gg = -0.001;//重力加速度

//速度矢量
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 lambda[] = { -5. * 4 / 12, 1. / 3, 1. / 12 };
//定义相关物理量
double rho[NX + 1][NY + 1]/*密度*/,
u[NX + 1][NY + 1][2]/*n+1时层的速度*/,
u0[NX + 1][NY + 1][2],//n时层的速度
t[NX + 1][NY + 1], //n+1时层的温度
p[NX + 1][NY + 1],//动压力
g[NX + 1][NY + 1][Q], //演化前的温度分布函数
G[NX + 1][NY + 1][Q], //演化后的温度分布函数
f[NX + 1][NY + 1][Q], //演化前的耦合分布函数
F[NX + 1][NY + 1][Q], //演化后的耦合分布函数
force; //外力项


int i, j, k, ip, jp, n;
double c, Re, Ra, Tm, DT, dx, dy, Lx, Ly, dt, rho0, Pr, Ma, niu, X, tau_f, tau_T, error, BETA, cs, uc;
double tempt, tempu, tempv;

//定义子函数
void init();
double feq(int k, double p, double u,double v);
double geq(int k, double T, double u,double v);
void evolution();
void output(int m);
void Error();


//定义主函数
int main()
{
	using namespace std;
	init();
	for (n = 0;; n++)
	{
		evolution();
		if (n % 10 == 0)
		{
			Error();
			cout << "The " << n << " th computation result:" << endl
				<< "The u,v,T of point (NX/2,NY/2) is:" << setprecision(6)
				<< u[NX / 2, NY / 2][0] << " , " << u[NX / 2][NY / 2][1] << " , " << t[NX / 2][NY / 2] << endl;
			cout << "The max relative error of uvT is:"
				<< resetiosflags(ios::scientific) << error << endl;
			if (n >= 20)
			{
				if (n % 20 == 0)
					output(n);
				if (error < 1.0e-8)
					break;
			}
		}
	}
	return 0;
}


//定义物理场
void init()
{
	dx = 1.0;
	dy = 1.0;
	Lx = dx*double(NX);
	Ly = dy*double(NY);
	dt = dx;
	c = dx / dt;
	rho0 = 1.17;//与Tm对应的空气的密度
	Pr = 0.71;
	DT = TH - TL;
	Tm = (TH + TL) / 2;
	Ra = 1.0e4;
	Ma = 0.1;
	cs = c / sqrt(3);
	uc = Ma*cs;
	BETA = uc*uc / (-Gg*DT*Ly);
	tau_f = 1.0 / 2 + (Ma*Ly*sqrt(3 * Pr)) / (c*dt*sqrt(Ra));
	niu = c*c*(tau_f - 1.0 / 2)*dt / 3;
	X = niu / Pr;
	tau_T = 0.5 + 2 * X / (c*c*dt);
	cout << "tau_f=" << tau_f
		<< "tau_T=" << tau_T << endl;

	for (i = 0; i <=NX; i++)
	for (j = 0; j <=NY; j++)
	{
		u[i][j][0] = 0;
		u[i][j][1] = 0;
		rho[i][j] = rho0;
		p[i][j] = 1;
		t[i][j] = Tm;
		t[0][j] = TH;
		t[NX][j] = TL;
		for (k = 0; k < Q; k++)
		{
			f[i][j][k] = feq(k, p[i][j], u[i][j][0],u[i][j][1]);
		}
		for (k = 1; k <= 4; k++)
		{
			g[i][j][k] = geq(k, t[i][j], u[i][j][0],u[i][j][1]);
		}
	}
}

//计算速度场平衡态分布函数
double feq(int k, double p, double u,double v)
{
	double eu, uv, feq;
	eu = (e[k][0] * u + e[k][1] *v);
	uv = u * u + u * u;
	if (k == 0)
	{
		feq=lambda[0] * p+w[k] * (3. * eu + 4.5*eu*eu - 1.5*uv);
	}
	else if (k==1||k==2||k==3||k==4)
	{
		feq=lambda[1] * p+w[k] * (3. * eu + 4.5*eu*eu - 1.5*uv);
	}
	else
	{
		feq=lambda[2] * p +w[k] * (3. * eu + 4.5*eu*eu - 1.5*uv);
	}
	return feq;
}

//计算温度场平衡态分布函数
double geq(int k, double T, double u,double v)
{
	double eu,geq;
	eu = e[k][0] * u + e[k][1] * v;
	geq = T / 4.*(1. + 2. * eu);
	return geq;
}


//演化
void evolution()
{
	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];
			if (k == 2 || k == 4)
			{
				force = -e[k][1] * Gg*BETA*(t[ip][jp] - Tm)/2.;
				F[i][j][k] = f[ip][jp][k] + (feq(k, p[ip][jp], u[ip][jp][0],u[ip][jp][1]) - f[ip][jp][k]) / tau_f + force*dt;
			}
			else
			{
				force = 0;
				F[i][j][k] = f[ip][jp][k] + (feq(k, p[ip][jp], u[ip][jp][0],u[ip][jp][1]) - f[ip][jp][k]) / tau_f + force*dt;
			}
		}

		for (k = 1; k <= 4; k++)
		{
			ip = i - e[k][0];
			jp = j - e[k][1];
			//温度分布函数演化方程
			G[i][j][k] = g[ip][jp][k] + (geq(k, t[ip][jp], u[ip][jp][0],u[ip][jp][1]) - g[ip][jp][k]) / tau_T;
		}
	}

	//====================边界处理-----采用非平衡外推格式=========================
	//左右边界
	for (j = 0; j <=NY; j++)
	{
		t[0][j] = TH;
		u[0][j][0] = 0;
		u[0][j][1] = 0;
		p[0][j] = p[1][j];

		t[NX][j] = TL;
		u[NX][j][0] = 0;
		u[NX][j][1] = 0;
		p[NX][j] = p[NX - 1][j];

		for (k = 0; k < Q; k++)
		{
			f[0][j][k] = feq(k, p[1][j], u[0][j][0],u[0][j][1]) + (f[1][j][k] - feq(k, p[1][j], u[1][j][0],u[1][j][1]));
			f[NX][j][k] = feq(k, p[NX - 1][j], u[NX][j][0],u[NX][j][1]) + (f[NX - 1][j][k] - feq(k, p[NX - 1][j], u[NX - 1][j][0],u[NX-1][j][1]));
		}
		for (k = 1; k <= 4; k++)
		{
			g[0][j][k] = geq(k, t[0][j], u[0][j][0],u[0][j][1]) + (g[1][j][k] - geq(k, t[1][j], u[1][j][0],u[1][j][1]));
			g[NX][j][k] = geq(k, t[NX][j], u[NX][j][0],u[NX][j][1]) + (g[NX - 1][j][k] - geq(k, t[NX - 1][j], u[NX - 1][j][0],u[NX-1][j][1]));
		}

	}

	//上下边界
	for (i = 0; i < NX; i++)
	{
		//下边界
		t[i][0] = t[i][1];
		u[i][0][0] = 0;
		u[i][0][1] = 0;
		p[i][0] = p[i][1];

		//上边界
		t[i][NY] = t[i][NY - 1];
		u[i][NY][0] = 0;
		u[i][NY][1] = 0;
		p[i][NY] = p[i][NY - 1];
		for (k = 0; k < Q; k++)
		{
			f[i][0][k] = feq(k, p[i][0], u[i][0][0],u[i][0][1]) + f[i][1][k] - feq(k, p[i][1], u[i][1][0],u[i][1][1]);
			f[i][NY][k] = feq(k, p[i][NY], u[i][NY][0],u[i][NY][1]) + f[i][NY - 1][k] - feq(k, p[i][NY - 1], u[i][NY - 1][0],u[i][NY-1][1]);
		}
		for (k = 1; k <= 4; k++)
		{
			g[i][0][k] = geq(k, t[i][0], u[i][0][0],u[i][0][1]) + g[i][1][k] - geq(k, t[i][1], u[i][1][0],u[i][1][1]);
			g[i][NY][k] = geq(k, t[i][NY], u[i][NY][0],u[i][NY][1]) + g[i][NY - 1][k] - geq(k, t[i][NY - 1], u[i][NY - 1][0],u[i][NY-1][1]);
		}
	}

	//四个角点
	f[0][0][5] = f[1][1][7];
	g[0][0][5] = g[1][1][7];
	f[NX][0][6] = f[NX - 1][1][8];
	g[NX][0][6] = g[NX - 1][1][8];
	f[0][NY][8] = f[1][NY - 1][6];
	g[0][NY][8] = g[1][NY - 1][6];
	f[NX][NY][7] = f[NX - 1][NY - 1][5];
	g[NX][NY][7] = g[NX - 1][NY - 1][5];

	//=======================================================================================

	//计算宏观量
	for (i = 1; i < NX; i++)
	for (j = 1; j < NY; j++)
	{
		u0[i][j][0] = u[i][j][0];
		u0[i][j][1] = u[i][j][1];
		p[i][j] = 0;
		tempt = 0;
		tempu = 0;
		tempv = 0;
		for (k = 0; k < Q; k++)
		{
			f[i][j][k] = F[i][j][k];
			g[i][j][k] = G[i][j][k];
			tempt += g[i][j][k];
			tempu += e[k][0] * f[i][j][k];
			tempv += e[k][1] * f[i][j][k];
		}
		t[i][j] = tempt;
		u[i][j][0] =tempu;
		u[i][j][1] =tempv;
		p[i][j] = (tempt - 1.5*(u[i][j][0] * u[i][j][0] + u[i][j][1] * u[i][j][1])) /(-lambda[0]);
	}
}


//输出整个流场的计算数据
void output(int m)
{
	ostringstream name;
	name << " Closed square cavity natural convection_" << m << ".plt";
	ofstream out(name.str().c_str());
	out << "Title=\"LBM  Closed square cavity natural convection\"\n" << "Variables=\"x\",\"y\",\"u\",\"v\",\"T\"\n"
		<< "ZONE T=\"BOX\",I="
		<< NX + 1 << ",J=" << NY + 1 << ",F=POINT" << endl;
	for (j = 0; j <= NY; j++)
	for (i = 0; i <= NX; i++)
	{
		out << double(i) / Lx << " "
			<< double(j) / Ly << " "
			<< u[i][j][0] << " " << u[i][j][1] << " " << t[i][j]  << endl;
	}
}


//计算相邻时层的温度和速度的最大相对误差
void Error()
{
	double temp1, temp2;
	temp1 = 0;
	temp2 = 0;
	for (i = 1; i < NX; i++)
	for (j = 1; j < NY; j++)
	{
		temp1 += (
			(u[i][j][0] - u0[i][j][0])*(u[i][j][0] - u0[i][j][0]) +
			(u[i][i][1] - u0[i][j][1])*(u[i][j][1] - u0[i][j][1]) 
			);
		temp2 += (
			u[i][j][0] * u[i][j][0] +
			u[i][j][1] * u[i][j][1] 
			);
	}
	temp1 = sqrt(temp1);
	temp2 = sqrt(temp2);
	error = temp1 / (temp2 + 1.0e-20);
}