d2q9: implement Zou/He boundaries on all sides

only the top wall can be imprinted with a velocity,
for the other sides the velocities are zero.

diff --git a/d2q9/esrc/d2q9.c b/d2q9/esrc/d2q9.c
index bc67bee41b20a728dfa3c1c1eac78530a5620108..bb7a940419fbc293623bccabff138fc4489ae64c 100644 (file)
--- a/d2q9/esrc/d2q9.c
+++ b/d2q9/esrc/d2q9.c
@@ -37,6 +37,44 @@ void d2q9_init(d2q9_block_t block)
 
 void d2q9_collide(d2q9_block_t f, int x, int y, FLOAT omega)
 {
+       if(row == 0 && y == 0) {
+               /* Zou/He boundary at top, with velocity */
+               const FLOAT UX = 0.0;           /* speed in x-direction */
+               const FLOAT UY = 0.0;           /* speed in y-direction */
+
+               FLOAT rho = ( f[y][x][0] + f[y][x][6] + f[y][x][2] +
+                       2 * ( f[y][x][5] + f[y][x][3] + f[y][x][4] ) ) /
+                       (1 - UY);
+
+               FLOAT tmp1 = ( f[y][x][6] - f[y][x][2] ) / 2;
+               FLOAT tmp2 = rho * (-UX - UY / 3) / 2;
+
+               f[y][x][7] = f[y][x][3] - tmp1 - tmp2;
+               f[y][x][1] = f[y][x][5] + tmp1 + tmp2;
+               f[y][x][8] = f[y][x][4] + 2 * rho * UY / 3;
+
+       } else if(row == CORES_Y-1 && y == BLOCKS_Y-1) {
+               /* Zou/He boundary at bottom, no velocity */
+               FLOAT tmp = ( f[y][x][6] - f[y][x][2] ) / 2;
+               f[y][x][3] = f[y][x][7] + tmp;
+               f[y][x][5] = f[y][x][1] - tmp;
+               f[y][x][4] = f[y][x][8];
+
+       } else if(col == 0 && x == 0) {
+               /* Zou/He boundary at left, no velocity */
+               FLOAT tmp = ( f[y][x][8] - f[y][x][4] ) / 2;
+               f[y][x][5] = f[y][x][1] + tmp;
+               f[y][x][7] = f[y][x][3] - tmp;
+               f[y][x][6] = f[y][x][2];
+
+       } else if(col == CORES_X-1 && x == BLOCKS_X-1) {
+               /* Zou/He boundary at right, no velocity */
+               FLOAT tmp = ( f[y][x][8] - f[y][x][4] ) / 2;
+               f[y][x][1] = f[y][x][5] - tmp;
+               f[y][x][3] = f[y][x][7] + tmp;
+               f[y][x][2] = f[y][x][6];
+       }
+
        /* macroscopic */
        FLOAT rho = f[y][x][0] + f[y][x][1] + f[y][x][2] +
                f[y][x][3] + f[y][x][4] + f[y][x][5] +