I am trying to simulate channel flow with D2Q9 LBM.

bbbbbbbb

i o

i o

i o

bbbbbbbb

where i is uniform velocity inlet

b -> wall boundary

o ->open boundary condition

I took this code to calculate stream function from one program, but not getting proper result.

strf(0,0)=0.0

do i=0,n

rhoav=0.5*(rho(i-1,0)+rho(i,0))

if(i.ne.0) strf(i,0)=strf(i-1,0)-rhoav*0.5*(v(i-1,0)+v(i,0))

do j=1,m

rhom=0.5*(rho(i,j)+rho(i,j-1))

strf(i,j)=strf(i,j-1)+rhom*0.5*(u(i,j-1)+u(i,j))

end do

end do

Please help me to figure out mistake.