PBCs.c.rar_PBCs

/* This is a standard two dimensional lattice Boltzmann program with periodic boundary conditions that should perform reasonably efficiently. It has been written for teaching purposes by Alexander Wagner at NDSU (March 2003). */ #include <stdio.h> #include <stdlib.h> #include <string.h> #include <unistd.h> /* for the sleep function. this may not be standard */ #include <math.h> #include <mygraph.h> #define xdim 31 #define ydim 21 #define DEBUG double f0[xdim][ydim], f1[xdim][ydim],f2[xdim][ydim],f3[xdim][ydim],f4[xdim][ydim], f5[xdim][ydim],f6[xdim][ydim],f7[xdim][ydim],f8[xdim][ydim]; #ifdef DEBUG double b1[xdim][ydim],b3[xdim][ydim],b5[xdim][ydim],b6[xdim][ydim], bb[xdim][ydim]; #endif double omega=1; /* Some constants that appear often and don't need to be calculated over and over again. So we calculate them just once here. */ double fourOnine=4.0/9.0, oneOnine=1.0/9.0, oneOthirtysix=1.0/36.0; /* Some variables for the suspended particle */ typedef double *doublep; typedef doublep pair[2]; pair *Bf1=NULL,*Bf3=NULL,*Bf5=NULL,*Bf6=NULL; int Bf1c,Bf1m=0,Bf3c,Bf3m=0,Bf5c,Bf5m=0,Bf6c,Bf6m=0; double Cx=15,Cy=10,Mx=0,My=0,Zxx=0,Zxy=0,Zyy=0,Ux=0,Uy=0,R=5,M=100; int Repeat=1,iterations=0,FrequencyMeasure=100,Graphics=1; int Pause=1,sstep=0,done=0; double Amp=0.001,Amp2=0.1,Width=10; /* Some special data types that are required to view the graphics */ int Xdim=xdim,Ydim=ydim,noreq; double rho[xdim][ydim]; int rhoreq=0; double u[xdim][ydim][2]; int ureq=0; void BCmem(pair **p, int c, int *cm){ if (c<*cm-2) return; *cm+=100; *p=(pair *) realloc(*p,*cm*sizeof(pair)); } void circleBC(){ /* Sets up the links that are cut by a circular object. */ int x,y,Px,Py,Pxp,Pyp,R2; double dx,dy; #ifdef DEBUG for (x=0;x<xdim;x++) for (y=0;y<ydim;y++) b1[x][y]=b3[x][y]=b5[x][y]=b6[x][y]=bb[x][y]=0; #endif Bf1c=Bf3c=Bf5c=Bf6c=0; R2=R*R; for (x=Cx-R-2;x<Cx+R+2;x++){ Px=(x+xdim)%xdim; Pxp=(x+1+xdim)%xdim; dx=x-Cx; for (y=Cy-R-2;y<Cy+R+2;y++){ Py=(y+ydim)%ydim; Pyp=(y+1+ydim)%ydim; dy=y-Cy; #ifdef DEBUG if ((dx*dx+dy*dy-R2)<=0) bb[Px][Py]=1;else bb[Px][Py]=-1; #endif if ((dx*dx+dy*dy-R2)*((dx+1)*(dx+1)+dy*dy-R2)<=0){ BCmem(&Bf1,Bf1c,&Bf1m); Bf1[Bf1c][0]=&(f2[Px][Py]); Bf1[Bf1c++][1]=&(f1[(Pxp)][Py]); #ifdef DEBUG b1[Px][Py]=1; } else { b1[Px][Py]=-1; #endif } if ((dx*dx+dy*dy-R2)*((dx)*(dx)+(dy+1)*(dy+1)-R2)<=0){ BCmem(&Bf3,Bf3c,&Bf3m); Bf3[Bf3c][0]=&(f4[Px][Py]); Bf3[Bf3c++][1]=&(f3[Px][Pyp]); #ifdef DEBUG b3[Px][Py]=1; } else { b3[Px][Py]=-1; #endif } if ((dx*dx+dy*dy-R2)*((dx+1)*(dx+1)+(dy+1)*(dy+1)-R2)<=0){ BCmem(&Bf5,Bf5c,&Bf5m); Bf5[Bf5c][0]=&(f7[Px][Py]); Bf5[Bf5c++][1]=&(f5[Pxp][Pyp]); #ifdef DEBUG b5[Px][Py]=1; } else { b5[Px][Py]=-1; #endif } if (((dx)*(dx)+(dy+1)*(dy+1)-R2)*((dx+1)*(dx+1)+(dy)*(dy)-R2)<=0){ BCmem(&Bf6,Bf6c,&Bf6m); Bf6[Bf6c][0]=&(f8[Pxp][Py]); Bf6[Bf6c++][1]=&(f6[Px][Pyp]); #ifdef DEBUG b6[Px][Py]=1; } else { b6[Px][Py]=-1; #endif } } } } void initParticle(){ Cx=10;Cy=10;Mx=0;My=0;Ux=0;Uy=0;R=5;M=100; } void init(){ int x,y; double n,ux,uy,uxx,uyy,uxy,usq; printf("initialize\n"); iterations=0; for (x=0;x<xdim;x++) for (y=0;y<ydim;y++){ /* here we can define the macroscopic properties of the initial condition. */ n=1+Amp2*exp(-(pow(x-xdim/2,2)+pow(y-ydim/2,2))/Width); ux=0/*.05*sin((y-ydim/2)*2*M_PI/xdim)*/; uy=0; /*n=1;ux=0;uy=0;*/ /* The following code initializes the f to be the local equilibrium values associated with the density and velocity defined above.*/ uxx=ux*ux; uyy=uy*uy; uxy=2*ux*uy; usq=uxx+uyy; f0[x][y]=fourOnine*n*(1-1.5*usq); f1[x][y]=oneOnine*n*(1+3*ux+4.5*uxx-1.5*usq); f2[x][y]=oneOnine*n*(1-3*ux+4.5*uxx-1.5*usq); f3[x][y]=oneOnine*n*(1+3*uy+4.5*uyy-1.5*usq); f4[x][y]=oneOnine*n*(1-3*uy+4.5*uyy-1.5*usq); f5[x][y]=oneOthirtysix*n*(1+3*(ux+uy)+4.5*(uxx+uxy+uyy)-1.5*usq); f6[x][y]=oneOthirtysix*n*(1+3*(-ux+uy)+4.5*(uxx-uxy+uyy)-1.5*usq); f7[x][y]=oneOthirtysix*n*(1+3*(-ux-uy)+4.5*(uxx+uxy+uyy)-1.5*usq); f8[x][y]=oneOthirtysix*n*(1+3*(ux-uy)+4.5*(uxx-uxy+uyy)-1.5*usq); } } void TotMomentum(){ int x,y; double Momx,Momy; Momx=Momy=0; for (x=0;x<xdim;x++) for (y=0;y<ydim;y++){ Momx+=f1[x][y]-f2[x][y]+f5[x][y]-f6[x][y]-f7[x][y]+f8[x][y]; Momy+=f3[x][y]-f4[x][y]+f5[x][y]+f6[x][y]-f7[x][y]-f8[x][y]; } printf("MomF = (%e,%e), MomP = (%e,%e), MomT = (%e,%e)\n", Momx,Momy,M*Ux,M*Uy,Momx+M*Ux,Momy+M*Uy); } void iterationColloid(){ #define alp 1 /* the degree of implicitness */ double zxx,zxy,zyy; /* The inverse Z tensor */ Cx+=Ux+xdim; Cx=fmod(Cx,xdim); Cy+=Uy+ydim; Cy=fmod(Cy,ydim); zxx=(M+alp*Zyy)/((M+alp*Zxx)*(M+alp*Zyy)+alp*alp*Zxy*Zxy); zxy= alp*Zxy /((M+alp*Zxx)*(M+alp*Zyy)+alp*alp*Zxy*Zxy); zyy=(M+alp*Zxx)/((M+alp*Zxx)*(M+alp*Zyy)+alp*alp*Zxy*Zxy); Ux+=zxx*(Mx-Zxx*Ux-Zxy*Uy)+zxy*(My-Zxy*Ux-Zyy*Uy); Uy+=zxy*(Mx-Zxx*Ux-Zxy*Uy)+zyy*(My-Zxy*Ux-Zyy*Uy); #undef alp } void iteration(){ int x,y,i; register double tmp1,tmp2; register double n,ux,uy,uxx,uyy,uxy,usq,Fx,Fy,Fxx,Fyy,Fxy,Fsq; double f1y[ydim],f2y[ydim],f5y[ydim],f6y[ydim],f7y[ydim],f8y[ydim]; double f3x[xdim],f4x[xdim],f5x[xdim],f6x[xdim],f7x[xdim],f8x[xdim]; /* first we perform the collision step */ for (x=0;x<xdim;x++) for (y=0;y<ydim;y++){ n=f0[x][y]+f1[x][y]+f2[x][y]+f3[x][y]+f4[x][y] +f5[x][y]+f6[x][y]+f7[x][y]+f8[x][y]; ux=f1[x][y]-f2[x][y]+f5[x][y]-f6[x][y]-f7[x][y]+f8[x][y]; uy=f3[x][y]-f4[x][y]+f5[x][y]+f6[x][y]-f7[x][y]-f8[x][y]; ux/=n; uy/=n; uxx=ux*ux; uyy=uy*uy; uxy=2*ux*uy; usq=uxx+uyy; /* We now need to implement any forcing terms. The term included here is just an example. You should not calculate a constant force term in the interation routine, but outside where it only gets calculated once.*/ Fx=Amp*sin(y*2*M_PI/ydim); 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; f0[x][y]+=omega*(fourOnine*n*(1-1.5*usq)-f0[x][y]) -fourOnine*1.5*Fsq; f1[x][y]+=omega*(oneOnine*n*(1+3*ux+4.5*uxx -1.5*usq)-f1[x][y]) +oneOnine*(3*Fx+4.5*Fxx-1.5*Fsq); f2[x][y]+=omega*(oneOnine*n*(1-3*ux+4.5*uxx -1.5*usq)-f2[x][y]) +oneOnine*(-3*Fx+4.5*Fxx-1.5*Fsq); f3[x][y]+=omega*(oneOnine*n*(1+3*uy+4.5*uyy -1.5*usq)-f3[x][y]) +oneOnine*(3*Fy+4.5*Fyy-1.5*Fsq); f4[x][y]+=omega*(oneOnine*n*(1-3*uy+4.5*uyy -1.5*usq)-f4[x][y]) +oneOnine*(-3*Fy+4.5*Fyy-1.5*Fsq); f5[x][y]+=omega*(oneOthirtysix*n*(1+3*(ux+uy)+4.5*(uxx+uxy+uyy) -1.5*usq)-f5[x][y]) +oneOthirtysix*(3*(Fx+Fy)+4.5*(Fxx+Fxy+Fyy)-1.5*Fsq); f6[x][y]+=omega*(oneOthirtysix*n*(1+3*(-ux+uy)+4.5*(uxx-uxy+uyy) -1.5*usq)-f6[x][y]) +oneOthirtysix*(3*(-Fx+Fy)+4.5*(Fxx-Fxy+Fyy)-1.5*Fsq); f7[x][y]+=omega*(oneOthirtysix*n*(1+3*(-ux-uy)+4.5*(uxx+uxy+uyy) -1.5*usq)-f7[x][y]) +oneOthirtysix*(3*(-Fx-Fy)+4.5*(Fxx+Fxy+Fyy)-1.5*Fsq); f8[x][y]+=omega*(oneOthirtysix*n*(1+3*(ux-uy)+4.5*(uxx-uxy+uyy) -1.5*usq)-f8[x][y]) +oneOthirtysix*(3*(Fx-Fy)+4.5*(Fxx-Fxy+Fyy)-1.5*Fsq); } /* now we need to move the densities according to their velocities we are using periodic boundary conditions */ /* since we are only using one lattice, we need to save some data on the boundaries so that it does not get overwritten */ for (y=0;y<ydim;y++){ f1y[y]=f1[xdim-1][y]; f2y[y]=f2[0][y]; f5y[y]=f5[xdim-1][y]; f6y[y]=f6[0][y]; f7y[y]=f7[0][y]; f8y[y]=f8[xdim-1][y]; } for (x=0;x<xdim;x++){ f3x[x]=f3[x][ydim-1]; f4x[x]=f4[x][0]; f5x[x]=f5[x][ydim-1]; f6x[x]=f6[x][ydim-1]; f7x[x]=f7[x][0]; f8x[x]=f8[x][0]; } /* Now we can move the densiti