#include "setting.h"

__global__ void Collision(int *flag, LBtype *pres, LBtype *vell, LBtype *f0, LBtype *forc)
{
	GLOBAL_INDEX()

	if(flag[k]==0)
	{
		LBtype P,Ux,Uy;
		LBtype M[9];
		LBtype g[9];

		LBtype Fx=forc[k*2], Fy=forc[k*2+1];

		for(int kk=0;kk<9;kk++)
			g[kk]=f0[k*9+kk];

		Ux=(g[1]+g[5]+g[8]-g[3]-g[6]-g[7]+0.5*Fx)/rho;
		Uy=(g[2]+g[5]+g[6]-g[4]-g[7]-g[8]+0.5*Fy)/rho;
		P =(g[0]+g[1]+g[2]+g[3]+g[4]+g[5]+g[6]+g[7]+g[8])/3.0;
		pressure[k]=P;

		M[0]=   g[0]  +g[1]  +g[2]  +g[3]  +g[4]  +g[5]  +g[6]  +g[7]  +g[8];
		M[1]=-4*g[0]  -g[1]  -g[2]  -g[3]  -g[4]+2*g[5]+2*g[6]+2*g[7]+2*g[8];
		M[2]= 4*g[0]-2*g[1]-2*g[2]-2*g[3]-2*g[4]  +g[5]  +g[6]  +g[7]  +g[8];
		M[3]=          g[1]         -g[3]         +g[5]  -g[6]  -g[7]  +g[8];
		M[4]=       -2*g[1]       +2*g[3]         +g[5]  -g[6]  -g[7]  +g[8];
		M[5]=                 g[2]         -g[4]  +g[5]  +g[6]  -g[7]  -g[8];
		M[6]=              -2*g[2]       +2*g[4]  +g[5]  +g[6]  -g[7]  -g[8];
		M[7]=          g[1]  -g[2]  +g[3]  -g[4];
		M[8]=                                      g[5]  -g[6]  +g[7]  -g[8];	

		M[0]=1.00*( 3*P	                   -M[0]);
		M[1]=1.20*(-6*P+3*rho*(Ux*Ux+Uy*Uy)-M[1])+(1-0.5*1.2)*6*(Ux*Fx+Uy*Fy);
		M[2]=1.20*( 3*P-3*rho*(Ux*Ux+Uy*Uy)-M[2])-(1-0.5*1.2)*6*(Ux*Fx+Uy*Fy);
		M[3]=1.00*( rho*Ux                 -M[3])+(1-0.5*0.0)*Fx;
		M[4]=1.20*(-rho*Ux                 -M[4])-(1-0.5*1.2)*Fx;
		M[5]=1.00*( rho*Uy                 -M[5])+(1-0.5*0.0)*Fy;
		M[6]=1.20*(-rho*Uy                 -M[6])-(1-0.5*1.2)*Fy;
		M[7]=  nu*(rho*(Ux*Ux-Uy*Uy)       -M[7])+(1-0.5*nu)*2*(Ux*Fx-Uy*Fy);
		M[8]=  nu*(rho*Ux*Uy               -M[8])+(1-0.5*nu)*(Ux*Fy+Uy*Fx);

		f0[k*9]  =g[0]+(  M[0]  -M[1]  +M[2])/9.0;
		f0[k*9+1]=g[1]+(4*M[0]  -M[1]-2*M[2]+6*M[3]-6*M[4]              +9*M[7])/36.0;
		f0[k*9+2]=g[2]+(4*M[0]  -M[1]-2*M[2]              +6*M[5]-6*M[6]-9*M[7])/36.0;
		f0[k*9+3]=g[3]+(4*M[0]  -M[1]-2*M[2]-6*M[3]+6*M[4]              +9*M[7])/36.0;
		f0[k*9+4]=g[4]+(4*M[0]  -M[1]-2*M[2]              -6*M[5]+6*M[6]-9*M[7])/36.0;
		f0[k*9+5]=g[5]+(4*M[0]+2*M[1]  +M[2]+6*M[3]+3*M[4]+6*M[5]+3*M[6]       +9*M[8])/36.0;
		f0[k*9+6]=g[6]+(4*M[0]+2*M[1]  +M[2]-6*M[3]-3*M[4]+6*M[5]+3*M[6]       -9*M[8])/36.0;
		f0[k*9+7]=g[7]+(4*M[0]+2*M[1]  +M[2]-6*M[3]-3*M[4]-6*M[5]-3*M[6]       +9*M[8])/36.0;
		f0[k*9+8]=g[8]+(4*M[0]+2*M[1]  +M[2]+6*M[3]+3*M[4]-6*M[5]-3*M[6]       -9*M[8])/36.0;
	}
}

__global__ void Streaming(LBtype *f0, LBtype *f1)
{
	GLOBAL_INDEX()
	int neighbor,nex,ney;
	int e[9][2]={{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};

	for(int kk=0;kk<9;kk++)
	{
		nex=(x+e[kk][0]+NX)%NX;
		ney=(y+e[kk][1]+NY)%NY;
		neighbor=ney*NX+nex;
		f1[neighbor*9+kk]=f0[k*9+kk];
	}
}

__global__ void BounceBack(int *flag, LBtype *f0)
{
	GLOBAL_INDEX()
	int neighbor,nex,ney;
	int e[9][2]={{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};
	int opp[9]={0,3,4,1,2,7,8,5,6};

	if(flag[k]==1)
		for(int kk=1;kk<9;kk++)
		{
			nex=(x+e[kk][0]+NX)%NX;
			ney=(y+e[kk][1]+NY)%NY;
			neighbor=ney*NX+nex;
			if(flag[neighbor]==0)
				f0[neighbor*9+kk]=f0[k*9+opp[kk]];
		}
}