collision.c

一个关于LB的经典程序。希望大家有好的东东就那出来

#include "global.h"
void collision()
{register int x,y;
register double rho,/* density */  vx,/* x-velocity */	vy,/* y-velocity */	square,tau_inv,	
f_eq0,f_eq1,f_eq2,f_eq3,f_eq4,f_eq5,f_eq6,f_eq7,f_eq8,/* equilibrium distr. */			dummy,product;		double 	**f0=f[0],		**f1=f[1],			**f2=f[2],
			**f3=f[3],
			**f4=f[4],
			**f5=f[5],
			**f6=f[6],
			**f7=f[7],
			**f8=f[8];
	/* multiplication is better than division in the loops to follow */
	tau_invvx*rho;
				f_eq2=rho*(1. + 3.0*vy + 4.5 *vy*vy - square);
				f_eq4=f_eq2-6.0*vy*rho;
	
				rho*=0.25;

				product=vx+vy;
				f_eq5=rho* (1. + 3.0*product + 4.5 *product*product - square);
				f_eq7=f_eq5-6.0*product*rho;		{
				/* introduce velocity boundary conditions for Driven Cavity
				by setting equilibrium distributions (see above) for rho, u_0, and v=0 */
				rho=f0[y][x]+f1[y][x]+f2[y][x]+f3[y][x]+f4[y][x]+f5[y][x]+f6[y][x]+f7[y][x]+f8[y][x];
				//rho=rho_0;
				vx=u_0;
				vy=0.0;
				square =1.5*(vx * vx +vy *vy );
				f0[y][x] =4./9.*rho*(1. - square);
	
				rho*=0.1111111111111111111111;
				f1[y][x]=rho*(1. + 3.0*vx + 4.5 *vx*vx - square);
				f3[y][x]=f1[y][x]-6.0*vx*rho;
				f2[y][x]=rho*(1. + 3.0*vy + 4.5 *vy*vy - square);
				f4[y][x]=f2[y][x]-6.0*vy*rho;
	
				rho*=0.25;

				product=vx+vy;
				f5[y][x]=rho* (1. + 3.0*product + 4.5 *product*product - square);
				f7[y][x]=f5[y][x]-6.0*product*rho;
				product=-vx+vy;
				f6[y][x]=rho* (1. + 3.0*product + 4.5 *product*product - square);
				f8[y][x]=f6[y][x]-6.0*product*rho;
			}
			else /* bounce back */
			{
				/* swap antiparallel distributions */
				dummy=f1[y][x];
				f1[y][x]=f3[y][x];
				f3[y][x]=dummy;
				dummy=f2[y][x];
				f2[y][x]=f4[y][x];
				f4[y][x]=dummy;
				dummy=f5[y][x];
				f5[y][x]=f7[y][x];
				f7[y][x]=dummy;
				dummy=f6[y][x];
				f6[y][x]=f8[y][x];
				f8[y][x]=dummy;
			}
		}
	}
}