function f = stream_2( Nx, Ny, f, f_temp ) % This function streams the distribution function for 2-D lattice Boltzmann % codes. The distribution f_temp is streamed to update the distribution f. % The code asumes standard D2Q9 directions: % 6 2 5 % 3 9 1 % 7 4 8 % Note the stationary direction is denoted by 9. % Boundary conditions are not included in this function, so separate % subroutines are needed to take care of the boundaries. i = 1; % left boundary i_p = (i+1); j = 1; % bottom left-hand corner j_p = (j+1); f(3, i, j) = f_temp(3, i_p, j); f(4, i, j) = f_temp(4, i, j_p); f(7, i, j) = f_temp(7, i_p, j_p); f(9, i, j) = f_temp(9, i, j); for j = 2:Ny-1 j_n = (j-1); j_p = (j+1); f(2, i, j) = f_temp(2, i, j_n); f(3, i, j) = f_temp(3, i_p, j); f(4, i, j) = f_temp(4, i, j_p); f(6, i, j) = f_temp(6, i_p, j_n); f(7, i, j) = f_temp(7, i_p, j_p); f(9, i, j) = f_temp(9, i, j); end j = Ny; % top left-hand corner j_n = (j-1); f(2, i, j) = f_temp(2, i, j_n); f(3, i, j) = f_temp(3, i_p, j); f(6, i, j) = f_temp(6, i_p, j_n); f(9, i, j) = f_temp(9, i, j); for i = 2:Nx-1 i_n = (i-1); i_p = (i+1); j = 1; % bottom boundary j_p = (j+1); f(1, i, j) = f_temp(1, i_n, j); f(3, i, j) = f_temp(3, i_p, j); f(4, i, j) = f_temp(4, i, j_p); f(7, i, j) = f_temp(7, i_p, j_p); f(8, i, j) = f_temp(8, i_n, j_p); f(9, i, j) = f_temp(9, i, j); for j = 2:Ny-1 j_n = (j-1); j_p = (j+1); f(1, i, j) = f_temp(1, i_n, j); f(2, i, j) = f_temp(2, i, j_n); f(3, i, j) = f_temp(3, i_p, j); f(4, i, j) = f_temp(4, i, j_p); f(5, i, j) = f_temp(5, i_n, j_n); f(6, i, j) = f_temp(6, i_p, j_n); f(7, i, j) = f_temp(7, i_p, j_p); f(8, i, j) = f_temp(8, i_n, j_p); f(9, i, j) = f_temp(9, i, j); end j = Ny; % top boundary j_n = (j-1); f(1, i, j) = f_temp(1, i_n, j); f(2, i, j) = f_temp(2, i, j_n); f(3, i, j) = f_temp(3, i_p, j); f(5, i, j) = f_temp(5, i_n, j_n); f(6, i, j) = f_temp(6, i_p, j_n); f(9, i, j) = f_temp(9, i, j); end i = Nx; % right boundary i_n = (i-1); j = 1; % bottom right-hand corner j_p = (j+1); f(1, i, j) = f_temp(1, i_n, j); f(4, i, j) = f_temp(4, i, j_p); f(8, i, j) = f_temp(8, i_n, j_p); f(9, i, j) = f_temp(9, i, j); for j = 2:Ny-1 j_n = (j-1); j_p = (j+1); f(1, i, j) = f_temp(1, i_n, j); f(2, i, j) = f_temp(2, i, j_n); f(4, i, j) = f_temp(4, i, j_p); f(5, i, j) = f_temp(5, i_n, j_n); f(8, i, j) = f_temp(8, i_n, j_p); f(9, i, j) = f_temp(9, i, j); end j = Ny; % top right-hand corner j_n = (j-1); f(1, i, j) = f_temp(1, i_n, j); f(2, i, j) = f_temp(2, i, j_n); f(5, i, j) = f_temp(5, i_n, j_n); f(9, i, j) = f_temp(9, i, j); end