rendered paste body#include<iostream>
using std::cout;
using std::endl;
#include <vector>
using std::vector;
#include <fstream>
using std::ofstream;
#include <iomanip>
using std::setw;
#include <sstream>
#include <string>
using std::stringstream;
#include<cmath>
#include "site.hpp"
/*Lattice Properties and constants*/
int L;
double nu, omega, rho;
int maxTime, tSaveStep;
void setConstants(){
L = 100; /*Dimentions of the box detector*/
rho = 1;
nu = 1.0; /*Kinematic Viscosity */
/*Relaxation Parameter
omega = 1.0 / (3.0*nu+1.0/2.0); */
omega = 1.0;
maxTime = 1000; /*number of iterations */
tSaveStep = 20; /* frequency of periodic saves to disk*/
}
/*This function is to make stuff look cleaner*/
inline double XYtoL(int x, int y){return x*L+y;}
/*Some "constants" which we will use. This will speed up calcs */
const double Four_Nine = 4.0/9.0;
const double One_Nine = 1.0/9.0;
const double One_Thirtysix = 1.0/36.0;
/*Decalre memory etc for the lattice */
vector<site> lattice;
/*and a temporary storage for lattice*/
/*vector<site> lattice_temp;*/
vector<site> lattice_temp(lattice);
/*Function to construct the lattice*/
void ConstructLattice(int L /*size of Lattice */){
for (int x=0; x<L; ++x)
for (int y=0; y<L; ++y)
lattice.push_back(site(x,y));
}
/*Function to initialize the lattice with the default BGK values*/
void InitLattice(){
/*Initial Velocity Parameters of the Lattice*/
double ux=0;
double uy=0;
/*calculate x-components and y-components of velocity */
double uxx=ux*ux;
double uyy=uy*uy;
double uxy=2*ux*uy;
double usq=uxx+uyy;
vector<site>::iterator LatticeSite;
for(LatticeSite=lattice.begin(); LatticeSite!=lattice.end(); LatticeSite++){
if(!LatticeSite->SolidNode){
/* The following code initializes the f to be the local equilibrium
values associated with the density and velocity defined above.*/
LatticeSite->LBPar[0] = Four_Nine*rho*(1-1.5*usq);
LatticeSite->LBPar[1]=One_Nine*rho*(1+3*uy+4.5*uyy-1.5*usq);
LatticeSite->LBPar[2]=One_Nine*rho*(1+3*ux+4.5*uxx-1.5*usq);
LatticeSite->LBPar[3]=One_Nine*rho*(1-3*uy+4.5*uyy-1.5*usq);
LatticeSite->LBPar[4]=One_Nine*rho*(1-3*ux+4.5*uxx-1.5*usq);
LatticeSite->LBPar[5]=One_Thirtysix*rho*(1+3*(ux+uy)+4.5*(uxx+uxy+uyy)-1.5*usq);
LatticeSite->LBPar[6]=One_Thirtysix*rho*(1+3*(-ux+uy)+4.5*(uxx-uxy+uyy)-1.5*usq);
LatticeSite->LBPar[7]=One_Thirtysix*rho*(1+3*(-ux-uy)+4.5*(uxx+uxy+uyy)-1.5*usq);
LatticeSite->LBPar[8]=One_Thirtysix*rho*(1+3*(ux-uy)+4.5*(uxx-uxy+uyy)-1.5*usq);
} else {
/*This is the initial setting for a solid node*/
LatticeSite->LBPar[0]=0;
LatticeSite->LBPar[1]=0;
LatticeSite->LBPar[2]=0;
LatticeSite->LBPar[3]=0;
LatticeSite->LBPar[4]=0;
LatticeSite->LBPar[5]=0;
LatticeSite->LBPar[6]=0;
LatticeSite->LBPar[7]=0;
LatticeSite->LBPar[8]=0;
}
}
}
/*Define Geometry of the detector*/
void Geometry(){
int iter;
for (iter=0; iter<=L-1; iter++){
/*Left wall - x=0, y loop from 0 to L-1)*/
lattice[XYtoL(0,iter)].SolidNode=true;
/*Right wall - x=L-1, y loop from 0 to L-1)*/
lattice[XYtoL((L-1),iter)].SolidNode=true;
/*top wall - y=L-1, x loop from 0 to L-1)*/
lattice[XYtoL(iter,(L-1))].SolidNode=true;
/*bottom wall - y=0, x loop from 0 to L-1)*/
lattice[XYtoL(iter,0)].SolidNode=true;
}
}
/***************************************************************/
/**************BEGIN LB SIMULATION******************************/
/*The collision operator */
void Collision(){
/* Aloocate memory for forcing terms and velocity terms */
register double n,ux,uy,uxx,uyy,uxy,usq,Fx,Fy,Fxx,Fyy,Fxy,Fsq;
/*assign temp storage for the LBPars*/
register double LBP0,LBP1,LBP2,LBP3,LBP4,LBP5,LBP6,LBP7,LBP8;
/****************************/
/*Collisions*/
/****************************/
vector<site>::iterator CollisionIterator;
for(CollisionIterator=lattice.begin(); CollisionIterator!=lattice.end(); CollisionIterator++)
{
if(!CollisionIterator->SolidNode){
LBP0 = CollisionIterator->LBPar[0]; LBP1 = CollisionIterator->LBPar[1]; LBP2 = CollisionIterator->LBPar[2]; LBP3 = CollisionIterator->LBPar[3]; LBP4 = CollisionIterator->LBPar[4];
LBP5 = CollisionIterator->LBPar[5]; LBP6 = CollisionIterator->LBPar[6]; LBP7 = CollisionIterator->LBPar[7]; LBP8 = CollisionIterator->LBPar[8];
/*Calculate the number density of the site*/
n=LBP0+LBP1+LBP2+LBP3+LBP4+LBP5+LBP6+LBP7+LBP8;
/*calculate velocity components*/
ux=LBP2-LBP4+LBP5+LBP6-LBP7-LBP8;
ux/=n;
uy=LBP1-LBP3+LBP5+LBP8-LBP6-LBP7;
uy/=n;
//cout<<n<<endl;
/*Make the velocity components into parts we need */
uxx=ux*ux;
uyy=uy*uy;
uxy=2*ux*uy;
usq=uxx+uyy;
/****************************/
/*Forcing Terms*/
Fx=0;/*0.002*sin((CollisionIterator->y)*2*3.14159/L);*/
Fy=0;
Fxx=2*n*Fx*ux;
Fyy=2*n*Fy*uy;
Fxy=2*n*(Fx*uy+Fy*ux);
Fsq=Fxx+Fyy;
Fx*=n;
Fy*=n;
/****************************/
/*Update the probability density functional */
CollisionIterator->LBPar[0]+= omega*(Four_Nine*n*(1-1.5*usq)-LBP0)-Four_Nine*1.5*Fsq;
CollisionIterator->LBPar[1]+=omega*(One_Nine*n*(1+3*uy+4.5*uyy -1.5*usq)-LBP1)+One_Nine*(3*Fy+4.5*Fyy-1.5*Fsq);
CollisionIterator->LBPar[2]+=omega*(One_Nine*n*(1+3*ux+4.5*uxx -1.5*usq)-LBP2)+One_Nine*(3*Fx+4.5*Fxx-1.5*Fsq);
CollisionIterator->LBPar[3]+=omega*(One_Nine*n*(1-3*uy+4.5*uyy -1.5*usq)-LBP3)+One_Nine*(-3*Fy+4.5*Fyy-1.5*Fsq);
CollisionIterator->LBPar[4]+=omega*(One_Nine*n*(1-3*ux+4.5*uxx -1.5*usq)-LBP4)+One_Nine*(-3*Fx+4.5*Fxx-1.5*Fsq);
CollisionIterator->LBPar[5]+=omega*(One_Thirtysix*n*(1+3*(-ux+uy)+4.5*(uxx-uxy+uyy)-1.5*usq)-LBP5)+One_Thirtysix*(3*(-Fx+Fy)+4.5*(Fxx-Fxy+Fyy)-1.5*Fsq);
CollisionIterator->LBPar[6]+=omega*(One_Thirtysix*n*(1+3*(ux+uy)+4.5*(uxx+uxy+uyy)-1.5*usq)-LBP6)+One_Thirtysix*(3*(Fx+Fy)+4.5*(Fxx+Fxy+Fyy)-1.5*Fsq);
CollisionIterator->LBPar[7]+=omega*(One_Thirtysix*n*(1+3*(ux-uy)+4.5*(uxx-uxy+uyy)-1.5*usq)-LBP7)+One_Thirtysix*(3*(Fx-Fy)+4.5*(Fxx-Fxy+Fyy)-1.5*Fsq);
CollisionIterator->LBPar[8]+=omega*(One_Thirtysix*n*(1+3*(-ux-uy)+4.5*(uxx+uxy+uyy)-1.5*usq)-LBP8)+One_Thirtysix*(3*(-Fx-Fy)+4.5*(Fxx+Fxy+Fyy)-1.5*Fsq);
}
}
}
/*The propagation operator - with no slip boundary conditions / geometry */
void Propagation(){
/*we first need a temporary storage to make the new lattice
We just copy the original one to save information about boundaries*/
//vector<site> lattice_temp(lattice);
//lattice_temp = lattice;
register double temp_boundary;// x_new, y_new;
double heating = 1.0;
double cooling = 0.0;
//cout<<lattice_temp[XYtoL((3+1),2)].x;
vector<site>::iterator PropagationIterator;
for(PropagationIterator=lattice.begin(); PropagationIterator!=lattice.end(); PropagationIterator++)
{
if(!PropagationIterator->SolidNode){
if (PropagationIterator->y == 1) { PropagationIterator->LBPar[1]=heating; }
if (PropagationIterator->y == L-1) { PropagationIterator->LBPar[3]=cooling;}
lattice_temp[XYtoL((PropagationIterator->x),(PropagationIterator->y))].LBPar[0] = PropagationIterator->LBPar[0];
lattice_temp[XYtoL((PropagationIterator->x),((PropagationIterator->y)+1))].LBPar[1] = PropagationIterator->LBPar[1];
lattice_temp[XYtoL(((PropagationIterator->x)+1),(PropagationIterator->y))].LBPar[2] = PropagationIterator->LBPar[2];
lattice_temp[XYtoL((PropagationIterator->x),((PropagationIterator->y)-1))].LBPar[3] = PropagationIterator->LBPar[3];
lattice_temp[XYtoL(((PropagationIterator->x)-1),(PropagationIterator->y))].LBPar[4] = PropagationIterator->LBPar[4];
lattice_temp[XYtoL(((PropagationIterator->x)-1),((PropagationIterator->y)+1))].LBPar[5] = PropagationIterator->LBPar[5];
lattice_temp[XYtoL(((PropagationIterator->x)+1),((PropagationIterator->y)+1))].LBPar[6] = PropagationIterator->LBPar[6];
lattice_temp[XYtoL(((PropagationIterator->x)+1),((PropagationIterator->y)-1))].LBPar[7] = PropagationIterator->LBPar[7];
lattice_temp[XYtoL(((PropagationIterator->x)-1),((PropagationIterator->y)-1))].LBPar[8] = PropagationIterator->LBPar[8];
} else {
/*Bounceback Boundary Conditions */
temp_boundary=PropagationIterator->LBPar[1]; PropagationIterator->LBPar[1]=PropagationIterator->LBPar[3]; PropagationIterator->LBPar[3]=temp_boundary;
temp_boundary=PropagationIterator->LBPar[2]; PropagationIterator->LBPar[2]=PropagationIterator->LBPar[4]; PropagationIterator->LBPar[4]=temp_boundary;
temp_boundary=PropagationIterator->LBPar[5]; PropagationIterator->LBPar[5]=PropagationIterator->LBPar[7]; PropagationIterator->LBPar[7]=temp_boundary;
temp_boundary=PropagationIterator->LBPar[6]; PropagationIterator->LBPar[6]=PropagationIterator->LBPar[8]; PropagationIterator->LBPar[8]=temp_boundary;
/*And then propagate the boundary!*/
/*heat the bottom wall and cool top*/
//if (PropagationIterator->y == 0) { PropagationIterator->LBPar[1]=heating;} //PropagationIterator->LBPar[5]=heating*sin(45); PropagationIterator->LBPar[6]=heating*sin(45);}
//if (PropagationIterator->y == L-1) { PropagationIterator->LBPar[1]=cooling; PropagationIterator->LBPar[5]=cooling*sin(45); PropagationIterator->LBPar[6]=cooling*sin(45);}
/*and propagate, but to prevent overrun, dont propagate 0 values!*/
if (!PropagationIterator->LBPar[1]==0) lattice_temp[XYtoL((PropagationIterator->x),((PropagationIterator->y)+1))].LBPar[1] = PropagationIterator->LBPar[1];
if (!PropagationIterator->LBPar[2]==0) lattice_temp[XYtoL(((PropagationIterator->x)+1),(PropagationIterator->y))].LBPar[2] = PropagationIterator->LBPar[2];
if (!PropagationIterator->LBPar[3]==0) lattice_temp[XYtoL((PropagationIterator->x),((PropagationIterator->y)-1))].LBPar[3] = PropagationIterator->LBPar[3];
if (!PropagationIterator->LBPar[4]==0) lattice_temp[XYtoL(((PropagationIterator->x)-1),(PropagationIterator->y))].LBPar[4] = PropagationIterator->LBPar[4];
if (!PropagationIterator->LBPar[5]==0) lattice_temp[XYtoL(((PropagationIterator->x)-1),((PropagationIterator->y)+1))].LBPar[5] = PropagationIterator->LBPar[5];
if (!PropagationIterator->LBPar[6]==0) lattice_temp[XYtoL(((PropagationIterator->x)+1),((PropagationIterator->y)+1))].LBPar[6] = PropagationIterator->LBPar[6];
if (!PropagationIterator->LBPar[7]==0) lattice_temp[XYtoL(((PropagationIterator->x)+1),((PropagationIterator->y)-1))].LBPar[7] = PropagationIterator->LBPar[7];
if (!PropagationIterator->LBPar[8]==0) lattice_temp[XYtoL(((PropagationIterator->x)-1),((PropagationIterator->y)-1))].LBPar[8] = PropagationIterator->LBPar[8];
}
}
lattice_temp.swap(lattice);
}
/*Function to count momentum */
void CountMomentum(float* px, float* py){
double pf_x=0;
double pf_y=0;
vector<site>::iterator LatticeSite;
for(LatticeSite=lattice.begin(); LatticeSite!=lattice.end(); LatticeSite++){
pf_x+=LatticeSite->LBPar[2]-LatticeSite->LBPar[4]+LatticeSite->LBPar[5]+LatticeSite->LBPar[6]-LatticeSite->LBPar[7]-LatticeSite->LBPar[8];
pf_y+=LatticeSite->LBPar[1]-LatticeSite->LBPar[3]+LatticeSite->LBPar[5]+LatticeSite->LBPar[8]-LatticeSite->LBPar[6]-LatticeSite->LBPar[7];
}
/*
px = Ux*pf_x;
py = Uy*pf_y;
*/
//cout<<pf_y;
}
ofstream gplotscript("movie.gp");
void DumpFField(int series){
stringstream ss;
ss << "FField" << series << ".dat";
//std::string s = ss.str();
ofstream outfile(ss.str().c_str());
double pf_x=0;
double pf_y=0;
double pf=0;
vector<site>::iterator LatticeSite;
for(LatticeSite=lattice.begin(); LatticeSite!=lattice.end(); LatticeSite++){
pf_x=LatticeSite->LBPar[2]-LatticeSite->LBPar[4]-LatticeSite->LBPar[5]+LatticeSite->LBPar[6]+LatticeSite->LBPar[7]-LatticeSite->LBPar[8];
pf_y=LatticeSite->LBPar[1]-LatticeSite->LBPar[3]+LatticeSite->LBPar[5]-LatticeSite->LBPar[8]+LatticeSite->LBPar[6]-LatticeSite->LBPar[7];
pf=sqrt(pf_x*pf_x+pf_y*pf_y);
//outfile<<LatticeSite->x<<setw(15)<<LatticeSite->y<<setw(15)<<pf_x<<setw(15)<<pf_y<<endl;
if (!(LatticeSite->x ==0 || LatticeSite->y ==0 || LatticeSite->x == L-1 ||LatticeSite->y ==L-1)) {
outfile<<LatticeSite->x<<setw(15)<<LatticeSite->y<<setw(15)<<pf<<endl;
}
}
gplotscript<<"splot \"FField"<<series<<".dat\" with points palette pointtype 6 pointsize .2"<<endl;
gplotscript<<"set zrange restore"<<endl;
outfile.close();
}
int main(){
setConstants();
ConstructLattice(L);
Geometry();
InitLattice();
lattice_temp=lattice;
gplotscript<<"set view map"<<endl;
gplotscript<<"set palette negative"<<endl;
int stepnum;
for (int time=0; time<=maxTime; time++){
Collision();
Propagation();
if (time%tSaveStep==0){
cout<<time<<endl;
stepnum = time/tSaveStep;
DumpFField(stepnum);
}
}
return 0;
}