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 = 50; /*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 = 100; /*number of iterations */ tSaveStep = 20; /* frequency of periodic saves to disk*/} /*This macro is to make stuff look cleaner*/#define XYtoL(x, y)(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;/*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; /*Left wall - x=0, y loop from 0 to L-1)*/ for (iter=0; iter<=L-1; iter++) lattice[XYtoL(0,iter)].SolidNode=true; /*Right wall - x=L-1, y loop from 0 to L-1)*/ for (iter=0; iter<=L-1; iter++) lattice[XYtoL((L-1),iter)].SolidNode=true; /*top wall - y=L-1, x loop from 0 to L-1)*/ for (iter=0; iter<=L-1; iter++) lattice[XYtoL(iter,(L-1))].SolidNode=true; /*bottom wall - y=0, x loop from 0 to L-1)*/ for (iter=0; iter<=L-1; iter++) 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; /****************************/ /*Collisions*/ /****************************/ vector<site>::iterator CollisionIterator; for(CollisionIterator=lattice.begin(); CollisionIterator!=lattice.end(); CollisionIterator++) { if(!CollisionIterator->SolidNode){ /*Calculate the number density of the site*/ n=CollisionIterator->LBPar[0]+CollisionIterator->LBPar[1]+CollisionIterator->LBPar[2]+CollisionIterator->LBPar[3]+CollisionIterator->LBPar[4]+CollisionIterator->LBPar[5]+CollisionIterator->LBPar[6]+CollisionIterator->LBPar[7]+CollisionIterator->LBPar[8]; /*calculate velocity components*/ ux=CollisionIterator->LBPar[2]-CollisionIterator->LBPar[4]+CollisionIterator->LBPar[5]+CollisionIterator->LBPar[6]-CollisionIterator->LBPar[7]-CollisionIterator->LBPar[8]; ux/=n; uy=CollisionIterator->LBPar[1]-CollisionIterator->LBPar[3]+CollisionIterator->LBPar[5]+CollisionIterator->LBPar[8]-CollisionIterator->LBPar[6]-CollisionIterator->LBPar[7]; uy/=n; cout<<(CollisionIterator->LBPar[0])*ux<<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)-CollisionIterator->LBPar[0])-Four_Nine*1.5*Fsq; CollisionIterator->LBPar[1]+=omega*(One_Nine*n*(1+3*uy+4.5*uyy -1.5*usq)-CollisionIterator->LBPar[2])+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)-CollisionIterator->LBPar[1])+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)-CollisionIterator->LBPar[4])+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)-CollisionIterator->LBPar[3])+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)-CollisionIterator->LBPar[5])+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)-CollisionIterator->LBPar[6])+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)-CollisionIterator->LBPar[7])+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)-CollisionIterator->LBPar[8])+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); register double temp_boundary, x_new, y_new; //cout<<lattice_temp[XYtoL((3+1),2)].x; vector<site>::iterator PropagationIterator; for(PropagationIterator=lattice.begin(); PropagationIterator!=lattice.end(); PropagationIterator++) { if(!PropagationIterator->SolidNode){ 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!*/ /*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.swap(lattice_temp); }/*Function to count momentum */void CountMomentum(){ 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; }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; 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]; outfile<<LatticeSite->x<<setw(15)<<LatticeSite->y<<setw(15)<<pf_x<<setw(15)<<pf_y<<endl; } outfile.close(); }int main(){ setConstants(); ConstructLattice(L); Geometry(); InitLattice(); 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;}