lbm.c (21684B)
1 #include <stdio.h> 2 #include <stdlib.h> // dynamic allocation 3 #include <math.h> 4 5 // Floating point precision 6 //typedef float Float; 7 typedef double Float; 8 9 // 3D vector 10 typedef struct { 11 Float x; 12 Float y; 13 Float z; 14 } Float3; 15 16 17 //// SIMULATION PARAMETERS 18 19 // Number of dimensions 20 const int n = 3; 21 22 // Grid dims 23 //const unsigned int nx = 3; 24 //const unsigned int ny = 6; 25 //const unsigned int nz = 3; 26 const unsigned int nx = 37; 27 const unsigned int ny = 37; 28 const unsigned int nz = 37; 29 30 // Grid cell width 31 const Float dx = 1.0; 32 33 // Number of flow vectors in each cell 34 const int m = 19; 35 36 // Time step length 37 //const double dt = 1.0; 38 const double dt = 1.0e-3; 39 //const double dt = 1.0e-6; 40 41 // Simulation end time 42 //const Float t_end = 1.5e-4; 43 const double t_end = 2.0; 44 //const double t_end = 1.0; 45 //const Float t_end = 10.1; 46 47 const double t_file = 0.01; 48 49 // Fluid dynamic viscosity 50 const Float nu = 8.9e-4; 51 52 // Gravitational acceleration 53 //const Float3 g = {0.0, 0.0, -10.0}; 54 const Float3 g = {0.0, 0.0, 0.0}; 55 56 // Initial cell fluid density (dimensionless) 57 const Float rho0 = 1.0; 58 59 // Inital cell fluid velocity (dimensionless) 60 const Float3 u0 = {0.0, 0.0, 0.0}; 61 62 // Courant criteria limit 63 const Float C_max = 1.0; 64 65 66 //// FUNCTION DEFINITIONS 67 68 Float3 MAKE_FLOAT3(Float x, Float y, Float z) 69 { 70 Float3 v; 71 v.x = x; v.y = y; v.z = z; 72 return v; 73 } 74 75 // Dot product of two Float3 vectors 76 Float dot(Float3 a, Float3 b) 77 { 78 return a.x*b.x + a.y*b.y + a.z*b.z; 79 } 80 81 // Viscosity parameter 82 Float tau(void) { 83 return (6.0*nu*dt/(dx*dx) + 1.0)/2.0; 84 } 85 86 // Get i-th value from cell x,y,z 87 unsigned int idx( 88 const unsigned int x, 89 const unsigned int y, 90 const unsigned int z) 91 { 92 return x + nx*y + nx*ny*z; 93 } 94 95 // Get i-th value from cell x,y,z 96 unsigned int idxi( 97 const unsigned int x, 98 const unsigned int y, 99 const unsigned int z, 100 const unsigned int i) 101 { 102 return x + ((y + z*ny)*nx) + (nx*ny*nz*i); 103 } 104 105 // Get i-th weight 106 Float w(unsigned int i) 107 { 108 if (n == 3 && m == 19) { 109 if (i == 0) 110 return 1.0/3.0; 111 else if (i > 0 && i < 7) 112 return 1.0/18.0; 113 else 114 return 1.0/36.0; 115 } else { 116 fprintf(stderr, "Error in w: m = %d != 19", m); 117 fprintf(stderr, ", n = %d != 3\n", n); 118 exit(EXIT_FAILURE); 119 } 120 } 121 122 void set_e_values(Float3 *e) 123 { 124 if (n == 3 && m == 19) { 125 e[0] = MAKE_FLOAT3( 0.0, 0.0, 0.0); // zero vel. 126 e[1] = MAKE_FLOAT3( 1.0, 0.0, 0.0); // face +x 127 e[2] = MAKE_FLOAT3(-1.0, 0.0, 0.0); // face -x 128 e[3] = MAKE_FLOAT3( 0.0, 1.0, 0.0); // face +y 129 e[4] = MAKE_FLOAT3( 0.0,-1.0, 0.0); // face -y 130 e[5] = MAKE_FLOAT3( 0.0, 0.0, 1.0); // face +z 131 e[6] = MAKE_FLOAT3( 0.0, 0.0,-1.0); // face -z 132 e[7] = MAKE_FLOAT3( 1.0, 1.0, 0.0); // edge +x,+y 133 e[8] = MAKE_FLOAT3(-1.0,-1.0, 0.0); // edge -x,-y 134 e[9] = MAKE_FLOAT3(-1.0, 1.0, 0.0); // edge -x,+y 135 e[10] = MAKE_FLOAT3( 1.0,-1.0, 0.0); // edge +x,-y 136 e[11] = MAKE_FLOAT3( 1.0, 0.0, 1.0); // edge +x,+z 137 e[12] = MAKE_FLOAT3(-1.0, 0.0,-1.0); // edge -x,-z 138 e[13] = MAKE_FLOAT3( 0.0, 1.0, 1.0); // edge +y,+z 139 e[14] = MAKE_FLOAT3( 0.0,-1.0,-1.0); // edge -y,-z 140 e[15] = MAKE_FLOAT3(-1.0, 0.0, 1.0); // edge -x,+z 141 e[16] = MAKE_FLOAT3( 1.0, 0.0,-1.0); // edge +x,-z 142 e[17] = MAKE_FLOAT3( 0.0,-1.0, 1.0); // edge -y,+z 143 e[18] = MAKE_FLOAT3( 0.0, 1.0,-1.0); // edge +y,-z 144 } else { 145 fprintf(stderr, "Error in set_e_values: m = %d != 19", m); 146 fprintf(stderr, ", n = %d != 3\n", n); 147 exit(EXIT_FAILURE); 148 } 149 } 150 151 // Equilibrium distribution along flow vector e 152 Float feq( 153 Float rho, 154 Float w, 155 Float3 e, 156 Float3 u) 157 { 158 Float c2 = dx/dt; 159 return rho*w * (1.0 + 3.0/c2*dot(e,u) 160 + 9.0/(2.0*c2*c2)*dot(e,u)*dot(e,u) 161 - 3.0/(2.0*c2)*dot(u,u)*dot(u,u)); 162 } 163 164 // Initialize cell densities, velocities, and flow vectors 165 void init_rho_v(Float* rho, Float3* u) 166 { 167 unsigned int x, y, z; 168 for (z=0; z<nz; z++) { 169 for (y=0; y<ny; y++) { 170 for (x=0; x<nx; x++) { 171 172 // Set velocity to u0 173 u[idx(x,y,z)] = MAKE_FLOAT3(u0.x, u0.y, u0.z); 174 175 // Set density to rho0 176 rho[idx(x,y,z)] = rho0; 177 } 178 } 179 } 180 } 181 182 void init_f(Float* f, Float* f_new, Float* rho, Float3* u, Float3* e) 183 { 184 unsigned int x, y, z, i; 185 Float f_val; 186 187 for (z=0; z<nz; z++) { 188 for (y=0; y<ny; y++) { 189 for (x=0; x<nx; x++) { 190 for (i=0; i<m; i++) { 191 192 // Set fluid flow vectors to v0 193 f_val = feq(rho[idx(x,y,z)], w(i), e[i], u[idx(x,y,z)]); 194 f[idxi(x,y,z,i)] = f_val; 195 f_new[idxi(x,y,z,i)] = f_val; 196 } 197 } 198 } 199 } 200 } 201 202 // Bhatnagar-Gross-Kroop approximation collision operator 203 Float bgk( 204 Float f, 205 Float tau, 206 Float rho, 207 Float w, 208 Float3 e, 209 Float3 u) 210 { 211 // Without gravitational drag 212 //return f - (f - feq(rho, w, e, u))/tau; 213 214 // With gravitational drag 215 Float f_ext; // External force along e 216 Float m_f = dx*dx*dx*rho; // Fluid mass 217 Float3 f_g = {m_f*g.x, m_f*g.y, m_f*g.z}; // Gravitational force 218 f_ext = dot(f_g, e); // Drag force along e 219 return f - (f - feq(rho, w, e, u))/tau 220 + (2.0*tau - 1)/(2.0*tau)*3.0/w*f_ext; 221 } 222 223 // Cell fluid density 224 Float find_rho( 225 const Float* f, 226 const unsigned int x, 227 const unsigned int y, 228 const unsigned int z) 229 { 230 int i; 231 Float rho = 0.0; 232 for (i=0; i<m; i++) 233 rho += f[idxi(x,y,z,i)]; 234 return rho; 235 } 236 237 // Cell fluid velocity 238 Float3 find_u( 239 const Float* f, 240 const Float rho, 241 const Float3* e, 242 const unsigned int x, 243 const unsigned int y, 244 const unsigned int z) 245 { 246 Float3 u = {0.0, 0.0, 0.0}; 247 Float f_i; 248 unsigned int i; 249 for (i=0; i<m; i++) { 250 f_i = f[idxi(x,y,z,i)]; 251 u.x += f_i*e[i].x/rho; 252 u.y += f_i*e[i].y/rho; 253 u.z += f_i*e[i].z/rho; 254 } 255 256 // Check the Courant-Frederichs-Lewy condition 257 if ((u.x*dt/dx + u.y*dt/dx + u.z*dt/dx) > C_max) { 258 fprintf(stderr, "Error, the Courant-Friderichs-Lewy condition is not "); 259 fprintf(stderr, "satisfied.\nTry one or more of the following:\n"); 260 fprintf(stderr, "- Decrease the timestep (dt)\n"); 261 fprintf(stderr, "- Increase the cell size (dx)\n"); 262 fprintf(stderr, "- Decrease the fluid viscosity (nu)\n"); 263 fprintf(stderr, "- Decrease the fluid density (rho)\n"); 264 exit(EXIT_FAILURE); 265 } 266 267 return u; 268 } 269 270 // Lattice-Boltzmann collision step. 271 // Fluid distributions are modified towards the cell equilibrium. 272 // Values are read from f, and written to f, rho, and u. 273 void collide( 274 Float* f, 275 Float* rho, 276 Float3* u, 277 const Float3* e) 278 { 279 unsigned int x, y, z, i; 280 Float rho_new; 281 Float3 u_new; 282 283 // For each cell 284 /*#pragma omp parallel for default(none) \ 285 private(x, y, z, rho_new, u_new, i) \ 286 firstprivate(e) \ 287 shared(f, rho, u) \ 288 schedule(dynamic)*/ 289 for (z=0; z<nz; z++) { 290 for (y=0; y<ny; y++) { 291 for (x=0; x<nx; x++) { 292 293 // Calculate macroscopic parameters 294 rho_new = find_rho(f, x, y, z); 295 u_new = find_u(f, rho_new, e, x, y, z); 296 297 // Store macroscopic parameters 298 int idx_ = idx(x,y,z); 299 rho[idx_] = rho_new; 300 u[idx_] = u_new; 301 302 // Find new f values by fluid particle collision 303 for (i=0; i<m; i++) { 304 int idxi_ = idxi(x,y,z,i); 305 f[idxi_] = bgk(f[idxi_], tau(), rho_new, 306 w(i), e[i], u_new); 307 } 308 } 309 } 310 } 311 } 312 313 // Lattice-Boltzmann streaming step. 314 // Propagate fluid flows to cell neighbors. 315 // Boundary condition: Bounce back 316 void stream(Float* f, Float* f_new) 317 { 318 // For each cell 319 unsigned int x, y, z; 320 for (z=0; z<nz; z++) { 321 for (y=0; y<ny; y++) { 322 for (x=0; x<nx; x++) { 323 324 // Face 0 325 f_new[idxi(x,y,z,0)] = fmax(0.0, f[idxi(x, y, z, 0)]); 326 327 // Face 1 (+x): Bounce back 328 if (x < nx-1) 329 f_new[idxi(x+1, y, z, 1)] 330 = fmax(0.0, f[idxi(x, y, z, 1)]); 331 else 332 f_new[idxi( x, y, z, 2)] 333 = fmax(0.0, f[idxi(x, y, z, 1)]); 334 335 // Face 2 (-x): Bounce back 336 if (x > 0) 337 f_new[idxi(x-1, y, z, 2)] 338 = fmax(0.0, f[idxi(x, y, z, 2)]); 339 else 340 f_new[idxi( x, y, z, 1)] 341 = fmax(0.0, f[idxi(x, y, z, 2)]); 342 343 // Face 3 (+y): Bounce back 344 if (y < ny-1) 345 f_new[idxi( x,y+1, z, 3)] 346 = fmax(0.0, f[idxi(x, y, z, 3)]); 347 else 348 f_new[idxi( x, y, z, 4)] 349 = fmax(0.0, f[idxi(x, y, z, 3)]); 350 351 // Face 4 (-y): Bounce back 352 if (y > 0) 353 f_new[idxi( x,y-1, z, 4)] 354 = fmax(0.0, f[idxi(x, y, z, 4)]); 355 else 356 f_new[idxi( x, y, z, 3)] 357 = fmax(0.0, f[idxi(x, y, z, 4)]); 358 359 // Face 5 (+z): Bounce back 360 if (z < nz-1) 361 f_new[idxi( x, y,z+1, 5)] 362 = fmax(0.0, f[idxi(x, y, z, 5)]); 363 else 364 f_new[idxi( x, y, z, 6)] 365 = fmax(0.0, f[idxi(x, y, z, 5)]); 366 367 // Face 6 (-z): Bounce back 368 if (z > 0) 369 f_new[idxi( x, y,z-1, 6)] 370 = fmax(0.0, f[idxi(x, y, z, 6)]); 371 else 372 f_new[idxi( x, y, z, 5)] 373 = fmax(0.0, f[idxi(x, y, z, 6)]); 374 375 376 // Edge 7 (+x,+y): Bounce back 377 if (x < nx-1 && y < ny-1) 378 f_new[idxi(x+1,y+1, z, 7)] 379 = fmax(0.0, f[idxi(x, y, z, 7)]); 380 else if (x < nx-1) 381 f_new[idxi(x+1, y, z, 10)] 382 = fmax(0.0, f[idxi(x, y, z, 7)]); 383 else if (y < ny-1) 384 f_new[idxi( x,y+1, z, 9)] 385 = fmax(0.0, f[idxi(x, y, z, 7)]); 386 else 387 f_new[idxi( x, y, z, 8)] 388 = fmax(0.0, f[idxi(x, y, z, 7)]); 389 390 // Edge 8 (-x,-y): Bounce back 391 if (x > 0 && y > 0) 392 f_new[idxi(x-1,y-1, z, 8)] 393 = fmax(0.0, f[idxi(x, y, z, 8)]); 394 else if (x > 0) 395 f_new[idxi(x-1, y, z, 9)] 396 = fmax(0.0, f[idxi(x, y, z, 8)]); 397 else if (y > 0) 398 f_new[idxi( x,y-1, z, 10)] 399 = fmax(0.0, f[idxi(x, y, z, 8)]); 400 else 401 f_new[idxi( x, y, z, 7)] 402 = fmax(0.0, f[idxi(x, y, z, 8)]); 403 404 // Edge 9 (-x,+y): Bounce back 405 if (x > 0 && y < ny-1) 406 f_new[idxi(x-1,y+1, z, 9)] 407 = fmax(0.0, f[idxi(x, y, z, 9)]); 408 else if (x > 0) 409 f_new[idxi(x-1, y, z, 8)] 410 = fmax(0.0, f[idxi(x, y, z, 9)]); 411 else if (y < ny-1) 412 f_new[idxi( x,y+1, z, 7)] 413 = fmax(0.0, f[idxi(x, y, z, 9)]); 414 else 415 f_new[idxi( x, y, z, 10)] 416 = fmax(0.0, f[idxi(x, y, z, 9)]); 417 418 // Edge 10 (+x,-y): Bounce back 419 if (x < nx-1 && y > 0) 420 f_new[idxi(x+1,y-1, z, 10)] 421 = fmax(0.0, f[idxi(x, y, z, 10)]); 422 else if (x < nx-1) 423 f_new[idxi(x+1, y, z, 7)] 424 = fmax(0.0, f[idxi(x, y, z, 10)]); 425 else if (y > 0) 426 f_new[idxi( x,y-1, z, 8)] 427 = fmax(0.0, f[idxi(x, y, z, 10)]); 428 else 429 f_new[idxi( x, y, z, 9)] 430 = fmax(0.0, f[idxi(x, y, z, 10)]); 431 432 // Edge 11 (+x,+z): Bounce back 433 if (x < nx-1 && z < nz-1) 434 f_new[idxi(x+1, y,z+1, 11)] 435 = fmax(0.0, f[idxi(x, y, z, 11)]); 436 else if (x < nx-1) 437 f_new[idxi(x+1, y, z, 16)] 438 = fmax(0.0, f[idxi(x, y, z, 11)]); 439 else if (z < nz-1) 440 f_new[idxi( x, y,z+1, 15)] 441 = fmax(0.0, f[idxi(x, y, z, 11)]); 442 else 443 f_new[idxi( x, y, z, 12)] 444 = fmax(0.0, f[idxi(x, y, z, 11)]); 445 446 // Edge 12 (-x,-z): Bounce back 447 if (x > 0 && z > 0) 448 f_new[idxi(x-1, y,z-1, 12)] 449 = fmax(0.0, f[idxi(x, y, z, 12)]); 450 else if (x > 0) 451 f_new[idxi(x-1, y, z, 15)] 452 = fmax(0.0, f[idxi(x, y, z, 12)]); 453 else if (z > 0) 454 f_new[idxi( x, y,z-1, 16)] 455 = fmax(0.0, f[idxi(x, y, z, 12)]); 456 else 457 f_new[idxi( x, y, z, 11)] 458 = fmax(0.0, f[idxi(x, y, z, 12)]); 459 460 // Edge 13 (+y,+z): Bounce back 461 if (y < ny-1 && z < nz-1) 462 f_new[idxi( x,y+1,z+1, 13)] 463 = fmax(0.0, f[idxi(x, y, z, 13)]); 464 else if (y < ny-1) 465 f_new[idxi( x,y+1, z, 18)] 466 = fmax(0.0, f[idxi(x, y, z, 13)]); 467 else if (z < nz-1) 468 f_new[idxi( x, y,z+1, 17)] 469 = fmax(0.0, f[idxi(x, y, z, 13)]); 470 else 471 f_new[idxi( x, y, z, 14)] 472 = fmax(0.0, f[idxi(x, y, z, 13)]); 473 474 // Edge 14 (-y,-z): Bounce back 475 if (y > 0 && z > 0) 476 f_new[idxi( x,y-1,z-1, 14)] 477 = fmax(0.0, f[idxi(x, y, z, 14)]); 478 else if (y > 0) 479 f_new[idxi( x,y-1, z, 17)] 480 = fmax(0.0, f[idxi(x, y, z, 14)]); 481 else if (z > 0) 482 f_new[idxi( x, y,z-1, 18)] 483 = fmax(0.0, f[idxi(x, y, z, 14)]); 484 else 485 f_new[idxi( x, y, z, 13)] 486 = fmax(0.0, f[idxi(x, y, z, 14)]); 487 488 // Edge 15 (-x,+z): Bounce back 489 if (x > 0 && z < nz-1) 490 f_new[idxi(x-1, y,z+1, 15)] 491 = fmax(0.0, f[idxi(x, y, z, 15)]); 492 else if (x > 0) 493 f_new[idxi(x-1, y, z, 12)] 494 = fmax(0.0, f[idxi(x, y, z, 15)]); 495 else if (z < nz-1) 496 f_new[idxi( x, y,z+1, 11)] 497 = fmax(0.0, f[idxi(x, y, z, 15)]); 498 else 499 f_new[idxi( x, y, z, 16)] 500 = fmax(0.0, f[idxi(x, y, z, 15)]); 501 502 // Edge 16 (+x,-z) 503 if (x < nx-1 && z > 0) 504 f_new[idxi(x+1, y,z-1, 16)] 505 = fmax(0.0, f[idxi(x, y, z, 16)]); 506 else if (x < nx-1) 507 f_new[idxi(x+1, y, z, 11)] 508 = fmax(0.0, f[idxi(x, y, z, 16)]); 509 else if (z > 0) 510 f_new[idxi( x, y,z-1, 12)] 511 = fmax(0.0, f[idxi(x, y, z, 16)]); 512 else 513 f_new[idxi( x, y, z, 15)] 514 = fmax(0.0, f[idxi(x, y, z, 16)]); 515 516 // Edge 17 (-y,+z) 517 if (y > 0 && z < nz-1) 518 f_new[idxi( x,y-1,z+1, 17)] 519 = fmax(0.0, f[idxi(x, y, z, 17)]); 520 else if (y > 0) 521 f_new[idxi( x,y-1, z, 14)] 522 = fmax(0.0, f[idxi(x, y, z, 17)]); 523 else if (z < nz-1) 524 f_new[idxi( x, y,z+1, 13)] 525 = fmax(0.0, f[idxi(x, y, z, 17)]); 526 else 527 f_new[idxi( x, y, z, 18)] 528 = fmax(0.0, f[idxi(x, y, z, 17)]); 529 530 // Edge 18 (+y,-z) 531 if (y < ny-1 && z > 0) 532 f_new[idxi( x,y+1,z-1, 18)] 533 = fmax(0.0, f[idxi(x, y, z, 18)]); 534 else if (y < ny-1) 535 f_new[idxi( x,y+1, z, 13)] 536 = fmax(0.0, f[idxi(x, y, z, 18)]); 537 else if (z > 0) 538 f_new[idxi( x, y,z-1, 14)] 539 = fmax(0.0, f[idxi(x, y, z, 18)]); 540 else 541 f_new[idxi( x, y, z, 17)] 542 = fmax(0.0, f[idxi(x, y, z, 18)]); 543 544 } 545 } 546 } 547 } 548 549 // Swap Float pointers 550 void swapFloats(Float* a, Float* b) 551 { 552 Float* tmp = a; 553 a = b; 554 b = tmp; 555 } 556 557 // Print density values to file stream (stdout, stderr, other file) 558 void print_rho(FILE* stream, Float* rho) 559 { 560 unsigned int x, y, z; 561 for (z=0; z<nz; z++) { 562 for (y=0; y<ny; y++) { 563 for (x=0; x<nx; x++) { 564 fprintf(stream, "%f\t", rho[idx(x,y,z)]); 565 } 566 fprintf(stream, "\n"); 567 } 568 fprintf(stream, "\n"); 569 } 570 } 571 572 // Print velocity values from y-plane to file stream 573 void print_rho_yplane(FILE* stream, Float* rho, unsigned int y) 574 { 575 unsigned int x, z; 576 for (z=0; z<nz; z++) { 577 for (x=0; x<nx; x++) { 578 fprintf(stream, "%f\t", rho[idx(x,y,z)]); 579 } 580 fprintf(stream, "\n"); 581 } 582 } 583 584 585 // Print velocity values to file stream (stdout, stderr, other file) 586 void print_u(FILE* stream, Float3* u) 587 { 588 unsigned int x, y, z; 589 for (z=0; z<nz; z++) { 590 for (y=0; y<ny; y++) { 591 for (x=0; x<nx; x++) { 592 fprintf(stream, "%.1ex%.1ex%.1e\t", 593 u[idx(x,y,z)].x, 594 u[idx(x,y,z)].y, 595 u[idx(x,y,z)].z); 596 } 597 fprintf(stream, "\n"); 598 } 599 fprintf(stream, "\n"); 600 } 601 } 602 603 // Print velocity values from y-plane to file stream 604 void print_u_yplane(FILE* stream, Float3* u, unsigned int y) 605 { 606 unsigned int x, z; 607 for (z=0; z<nz; z++) { 608 for (x=0; x<nx; x++) { 609 fprintf(stream, "%.1ex%.1ex%.1e\t", 610 u[idx(x,y,z)].x, 611 u[idx(x,y,z)].y, 612 u[idx(x,y,z)].z); 613 } 614 fprintf(stream, "\n"); 615 } 616 } 617 618 619 int main(int argc, char** argv) 620 { 621 printf("### Lattice-Boltzman D%dQ%d test ###\n", n, m); 622 623 FILE* frho; 624 char filename[40]; 625 626 // Print parameter vals 627 //printf("Grid dims: nx = %d, ny = %d, nz = %d: %d cells\n", 628 //nx, ny, nz, ncells); 629 630 // Set cell flow vector values 631 Float3 e[m]; set_e_values(e); 632 633 // Particle distributions 634 unsigned int ncells = nx*ny*nz; 635 Float* f = malloc(ncells*m*sizeof(Float)); 636 Float* f_new = malloc(ncells*m*sizeof(Float)); 637 638 // Cell densities 639 Float* rho = malloc(ncells*sizeof(Float)); 640 641 // Cell flow velocities 642 Float3* u = malloc(ncells*sizeof(Float3)); 643 644 // Set densities, velocities and flow vectors 645 init_rho_v(rho, u); 646 rho[idx(nx/2,ny/2,nz/2)] *= 1.0001; 647 init_f(f, f_new, rho, u, e); 648 649 // Temporal loop 650 double t = 0.0; 651 double t_file_elapsed = 0.0; 652 653 // Save initial state 654 sprintf(filename, "out/rho_y%d_t%.2f.txt", ny/2, t); 655 if ((frho = fopen(filename, "w"))) { 656 print_rho_yplane(frho, rho, ny/2); 657 fclose(frho); 658 } else { 659 fprintf(stderr, "Error: Could not open output file "); 660 fprintf(stderr, "%s", filename); 661 fprintf(stderr, "\n"); 662 exit(EXIT_FAILURE); 663 } 664 665 // Temporal loop 666 for (t = 0.0; t < t_end; t += dt, t_file_elapsed += dt) { 667 668 // Report time to stdout 669 printf("\rt = %.1fs./%.1fs., %.1f%% done", t, t_end, t/t_end*100.0); 670 671 // LBM collision and streaming 672 collide(f, rho, u, e); 673 stream(f, f_new); 674 675 // Swap f and f_new 676 Float* tmp = f; 677 f = f_new; 678 f_new = tmp; 679 680 // Print x-z plane to file 681 if (t_file_elapsed >= t_file) { 682 sprintf(filename, "out/rho_y%d_t%.2f.txt", ny/2, t); 683 if ((frho = fopen(filename, "w"))) { 684 print_rho_yplane(frho, rho, ny/2); 685 fclose(frho); 686 } else { 687 fprintf(stderr, "Error: Could not open output file "); 688 fprintf(stderr, "%s", filename); 689 fprintf(stderr, "\n"); 690 exit(EXIT_FAILURE); 691 } 692 t_file_elapsed = 0.0; 693 } 694 } 695 printf("\n"); 696 697 // Report values to stdout 698 //fprintf(stdout, "rho\n"); 699 //print_rho(stdout, rho); 700 //fprintf(stdout, "u\n"); 701 //print_u(stdout, u); 702 703 // Clear memory 704 free(f); 705 free(f_new); 706 free(rho); 707 free(u); 708 709 return EXIT_SUCCESS; 710 }