Introduction
Introduction Statistics Contact Development Disclaimer Help
tlbm.c - lbm-d3q19 - 3D lattice-Boltzmann code to approximate Navier-Stokes inc…
git clone git://src.adamsgaard.dk/lbm-d3q19
Log
Files
Refs
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…
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…
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 }
You are viewing proxied material from mx1.adamsgaard.dk. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.