/* 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).
Extended to simulate the motion of one particle (March 2005).
*/
#include
#include
#include
#include /* for the sleep function. this may not be standard */
#include
#include
#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