tlbm.c - lbm-d3q19 - 3D lattice-Boltzmann code to approximate Navier-Stokes incompressible flow
HTML git clone git://src.adamsgaard.dk/lbm-d3q19
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tlbm.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 }