#ifdef __APPLE__ #include #else #include #endif #include #include #include #include #define WINDOW_WIDTH 100 #define WINDOW_HEIGHT 100 #define LATTICE_WIDTH 80 #define LATTICE_HEIGHT 80 #define VECTOR_SPACING 2 #define POINT_SIZE 4 #define SPACING 1 #define Q 9 #define DIM 2 #define EPS 0.1 double lattice[LATTICE_HEIGHT][LATTICE_WIDTH][Q]; double tmp_lattice[LATTICE_HEIGHT][LATTICE_WIDTH][Q]; double masses[LATTICE_HEIGHT][LATTICE_WIDTH]; double densities[LATTICE_HEIGHT][LATTICE_WIDTH]; int e[Q][DIM] = {{0,0}, {1, 0}, {0, 1}, {-1, 0}, {0, -1}, {1, 1}, {-1,1}, {-1,-1}, {1, -1}}; double weights[Q] = {4.0/9.0, 1.0/9.0 , 1.0/9.0 , 1.0/9.0 , 1.0/9.0 , 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0}; double viscocity = 9; double gravity = 9.8; double omega; int step = 0; double zeta = 4.0/5.0; double dt = 0; double get_omega(double viscocity) { return(2.0 / (6 * viscocity + 1)); } double get_rho(int x, int y) { int i; int r = 0; for(i = 0; i < Q; i++) { r += lattice[y][x][i]; } return r; } void get_u(int x, int y, double * u) { int i, j; for(i = 0; i < DIM; i++) { u[i] = 0; } for(i = 0; i < DIM; i++) { for(j = 0; j < Q; j++) { u[i] += lattice[y][x][j] * e[j][i]; } } } double get_u_squared(double * u) { int i; double u_squared = 0; for(i = 0; i < DIM; i++) { u_squared += u[i] * u[i]; } return u_squared; } double get_u_dot_e(int i, double * u) { int j; double u_dot_e_i = 0; for(j = 0; j < DIM; j++) { u_dot_e_i += e[i][j] * u[j]; } return u_dot_e_i; } double get_f_eq(int i, double rho, double * u) { double u_dot_e_i = get_u_dot_e(i, u); double u_squared = get_u_squared(u); return(weights[i] * (rho + 3 * u_dot_e_i - 1.5 * u_squared + 4.5 * u_dot_e_i * u_dot_e_i)); } double get_u_max(double * u) { int i; double u_max = 0; for(i = 0; i < Q; i++) { if(u_max < abs(u[i])) { u_max = abs(u[i]); } } return u_max; } double get_volume( ) { int x, y; double v = 0; for(y = 0; y < LATTICE_HEIGHT; y++) { for(x = 0; x< LATTICE_WIDTH; x++) { v += densities[y][x]; } } return v; } double get_mass( ) { int x, y; double m = 0; for(y = 0; y < LATTICE_HEIGHT; y++) { for(x = 0; x< LATTICE_WIDTH; x++) { m += masses[y][x]; } } return m; } double get_mass_difference(int x, int y) { int i; double mass = 0; for(i = 0; i < Q; i++ ) { mass += tmp_lattice[y][x][i] - lattice[y][x][i]; } return mass; } int resolve_collisions( ) { int x,y, i; int rho, rho_med; double u[DIM]; double f_eq; double u_max; for(x = 0; x < LATTICE_WIDTH; x++) { for(y = 0; y < LATTICE_HEIGHT; y++) { rho = get_rho(x, y); get_u(x, y, u); u_max = get_u_max(u); if(u_max > 1.0/(zeta * 6.0)) { dt *= zeta; double viscocity_old = viscocity; viscocity = zeta * (viscocity - 0.5) + 0.5; omega = get_omega(viscocity); gravity *= zeta * zeta; rho_med = get_volume( ) / get_mass( ); double rho_n, u_n[DIM], m_n, e_n; double s_f, s_t; for(x = 0; x < LATTICE_WIDTH; x++) { for(y = 0; y < LATTICE_HEIGHT; y++) { rho = get_rho(x, y); rho_n = zeta * (rho - rho_med) + rho_med; get_u(x,y, u); for(i = 0; i < DIM; i++) { u_n[i] = zeta * u[i]; } masses[y][x] *= rho / rho_n; densities[y][x] *= masses[y][x] / rho_n; for(i = 0; i < Q; i++) { s_f = get_f_eq(i, rho_n, u_n) / get_f_eq(i, rho, u); s_t = zeta * (viscocity / viscocity_old); tmp_lattice[y][x][i] = s_f * (get_f_eq(i, rho, u) + s_t * (lattice[y][x][i] - get_f_eq(i, rho,u))); } } } x = LATTICE_WIDTH; y = LATTICE_HEIGHT; } else { densities[y][x] = masses[y][x] / rho; for(i = 0; i < Q; i++) { f_eq = get_f_eq(i, rho, u); tmp_lattice[y][x][i] = (1.0 - omega) * lattice[y][x][i] + omega * f_eq; } masses[y][x] -= get_mass_difference(x,y); } } } return 1; } void propagate( ) { int x, y, i; int e_x, e_y; int new_x, new_y; for(x = 0; x < LATTICE_WIDTH; x++) { for(y = 0; y < LATTICE_HEIGHT; y++) { for(i = 0; i < Q; i++) { e_x = e[i][0]; e_y = e[i][1]; new_x = (x + e_x + LATTICE_WIDTH) % LATTICE_WIDTH; new_y = (y + e_y + LATTICE_HEIGHT) % LATTICE_HEIGHT; lattice[new_y][new_x][i] = tmp_lattice[y][x][i]; } } } } void next( ) { int run_propagate = resolve_collisions( ); if(run_propagate) { propagate( ); } } void display( ) { int x, y, i, j; int v_x, v_y; double slope[DIM]; glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); glBegin(GL_LINES); for (x = 0; x < LATTICE_WIDTH; x += VECTOR_SPACING) { for (y = 0; y < LATTICE_HEIGHT; y += VECTOR_SPACING) { for(j = 0; j < DIM; j++) { slope[j] = 0.0; } for(v_x = x; v_x < x + VECTOR_SPACING; v_x++) { for(v_y = y; v_y < y + VECTOR_SPACING; v_y++) { for(i = 0; i < Q; i++) { for(j = 0; j < DIM; j++) { slope[j] += e[i][j] * lattice[v_y][v_x][i]; } } } } for(j = 0; j < DIM; j++) { slope[j] /= 0.5 * (VECTOR_SPACING); } if(slope[DIM-1] < 0) { glColor3f(1.0, 0.0, 0.0); } else { glColor3f(0.0, 0.0, 1.0); } glVertex3f(x + WINDOW_WIDTH/8, y + WINDOW_HEIGHT/8, 0); glVertex3f(x + slope[0] + WINDOW_WIDTH/8, y + slope[1] + WINDOW_HEIGHT/8, 0); } } glEnd( ); glutSwapBuffers( ); } void idle( ) { if(step || !step) { next( ); step = 0; } glutPostRedisplay(); } void init_vars( ) { int i, j, k; omega = get_omega(viscocity); for (i = 0; i < LATTICE_WIDTH; i++) { for (j = 0; j < LATTICE_HEIGHT; j++) { for(k = 0; k < Q; k++) { lattice[i][j][k] = 0; tmp_lattice[i][j][k] = 0; masses[i][j] = 0; densities[i][j] = 0; } } } for (i = LATTICE_HEIGHT/3; i < 3*LATTICE_WIDTH/4; i++) { for (j = LATTICE_WIDTH/3; j < 3*LATTICE_HEIGHT/4; j++) { for(k = 0; k < Q; k++) { lattice[i][j][k] = 0.2; tmp_lattice[i][j][k] = 0.1; masses[i][j] = 0.1; densities[i][j] = 0.1; } } } } void init_gl() { glClearColor(0.0, 0.0, 0.0, 0.0); glMatrixMode(GL_PROJECTION); glLoadIdentity(); glOrtho(0.0, WINDOW_WIDTH, 0.0, WINDOW_HEIGHT, -20.0, 20.0); glMatrixMode(GL_MODELVIEW); glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST); glShadeModel(GL_SMOOTH); glPointSize(POINT_SIZE); glLineWidth(POINT_SIZE); } void keyboard(unsigned char key, int x, int y) { if(key == 'n') { step = 1; } } int main(int argc, char** argv) { glutInit(&argc, argv); glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH); glutInitWindowSize(WINDOW_WIDTH, WINDOW_HEIGHT); glutCreateWindow("D2Q9"); glEnable(GL_DEPTH_TEST); init_vars(); init_gl(); glutDisplayFunc(display); glutKeyboardUpFunc(keyboard); glutIdleFunc(idle); glutMainLoop(); return 0; }