www.pudn.com > cylindar.zip > cylindar.f90, change:2013-11-19,size:21061b


 
MODULE cellConst 
    INTEGER, PARAMETER:: fluid = 0, wall = 1, inlet = 10, outlet = 11, incylinder=100 
END MODULE cellConst 
 
 
!======================================================== 
! Lattice constants for the D2Q9 lattice 
!======================================================== 
MODULE D2Q9Const 
!	 D2Q9 Weights 
    REAL(8),PARAMETER:: t(0:8) = (/4.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0& 
                                           &,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0/) 
!	D2Q9 Directions 
    INTEGER:: v(0:8,0:1) 
!       = (/(/0,1,0,-1,0,1,-1,-1,1/),(/0,0,1,0,-1,1,1,-1,-1/)/) 
 
    INTEGER, PARAMETER:: opposite(0:8) = (/0,3,4,1,2,7,8,5,6/) 
END MODULE D2Q9Const 
 
 
!======================================================== 
! Constants for simulation setup 
!======================================================== 
MODULE simParam 
    !INTEGER, PARAMETER:: xdim = 4500 
    !INTEGER, PARAMETER:: ydim = 800 
    !INTEGER, PARAMETER:: obstX = xdim/5 
    !INTEGER, PARAMETER:: obstY = ydim/2 
    !INTEGER, PARAMETER:: obstR = ydim/10+1 
    INTEGER, PARAMETER:: tmax = 100000 
    REAL(8), PARAMETER:: umax = 0.2d0 
    REAL(8), PARAMETER:: Re1 = 500.0d0 
    INTEGER, PARAMETER:: xdim = 2200 
    INTEGER, PARAMETER:: ydim = 320 
    INTEGER, PARAMETER:: obstX = 150 
    INTEGER, PARAMETER:: obstY = 155 
    INTEGER, PARAMETER:: obstR = 50 
END MODULE simParam 
 
 
!======================================================== 
! The main program, implementing a flow past a cylinder 
!======================================================== 
 
PROGRAM unsteady 
    USE simParam, ONLY: xdim, ydim, tmax 
    USE D2Q9Const, ONLY: v 
    implicit none 
 
    REAL(8):: omega, time1, time2, timeTot 
    REAL(8), dimension(:,:,:), allocatable:: f, fEq, u, meq 
    REAL(8), dimension(:,:), allocatable:: rho, uSqr, streamline, vorticity 
    INTEGER, dimension(:,:), allocatable:: image 
    INTEGER:: tStep=0 
 
    allocate(f(ydim,xdim,0:8)) 
    allocate(fEq(ydim,xdim,0:8)) 
    allocate(meq(ydim,xdim,0:8)) 
    allocate(u(ydim,xdim,0:1)) 
    allocate(uSqr(ydim,xdim)) 
    allocate(rho(ydim,xdim)) 
    allocate(image(ydim,xdim)) 
    allocate(streamline(ydim,xdim)) 
    allocate(vorticity(ydim,xdim)) 
 
 
 
    CALL constructImage(image) 
    CALL computeOmega(omega) 
    CALL writeInput(omega) 
    CALL initMacro(rho,u,uSqr,image) 
    CALL computeFeq(fEq,rho,u,uSqr) 
    f = fEq 
    CALL computevorticity(u, vorticity, image) 
    !CALL computestreamline(u,streamline,image) 
    CALL writeOutput(u,tstep,vorticity) 
    timeTot = 0.0d0 
    do tStep = 1, tmax 
        CALL CPU_TIME(time1) 
        CALL inletOutlet(f,rho,u,image) 
        CALL boundaries(f,image) 
        CALL computeMacros(f,rho,u,uSqr) 
       ! CALL computeFeq(fEq,rho,u,uSqr) 
       ! CALL collide(f,fEq,omega,image) 
        CALL computemeq(meq,rho,u,uSqr) 
        CALL coll_mrt(f,meq,omega,image) 
        CALL streaming(f,v,image) 
!        CALL writeImage(image) 
        IF (MOD(tstep, 1000)==0) THEN 
            !CALL computestreamline(u,streamline,image) 
            CALL computevorticity(u, vorticity, image) 
            CALL writeOutput(u,tstep,vorticity) 
        END IF 
        CALL CPU_TIME(time2) 
        timeTot = timeTot + (time2-time1) 
    end do 
 
    !CALL computestreamline(u,streamline,image) 
    CALL computevorticity(u, vorticity, image) 
    CALL writeOutput(u,tmax,vorticity) 
 
    write(*,*) dble(tmax) * (dble(ydim * xdim)) / timeTot ,'cells per second' 
 
    deallocate(f) 
    deallocate(fEq) 
    deallocate(meq) 
    deallocate(u) 
    deallocate(uSqr) 
    deallocate(rho) 
    deallocate(image) 
    deallocate(streamline) 
    deallocate(vorticity) 
END PROGRAM unsteady 
 
 
!======================================================== 
! Compute the relaxation PARAMETER from the Reynolds number 
!======================================================== 
SUBROUTINE computeOmega(omega) 
    USE simParam, ONLY: Re1,umax,obstR 
 
    implicit none 
 
    REAL(8), INTENT(INOUT):: omega 
    REAL(8):: nu 
 
    nu    =  umax * 2.0d0 * dble(obstR) / Re1 
    omega = 1.0d0 / (3.0d0*nu+0.5d0) 
END SUBROUTINE computeOmega 
 
 
!======================================================== 
! Construct an array the defines the flow geometry 
!======================================================== 
SUBROUTINE constructImage(image) 
    USE cellConst 
    USE simParam, ONLY: xdim, ydim, obstX, obstY, obstR 
    USE D2Q9Const, ONLY: v 
 
    implicit none 
 
    INTEGER, INTENT(INOUT):: image(ydim,xdim) 
    INTEGER:: x,y 
 
    v(0:8,0) = (/0,1,0,-1,0,1,-1,-1,1/) 
    v(0:8,1) = (/0,0,1,0,-1,1,1,-1,-1/) 
 
    image          = fluid 
    image(:,1)     = inlet 
    image(:,xdim)  = outlet 
    image(1,:)     = wall 
    image(ydim,:)  = wall 
    do x=1,xdim 
      do y=1,ydim 
        if(((x-obstX)**2+(y-obstY)**2) <= (obstR**2)) image(y,x)=incylinder 
      enddo 
    enddo 
 
END SUBROUTINE constructImage 
 
 
!======================================================== 
! Initialize the simulation to Poiseuille profile at 
!       an equilibrium distribution 
!======================================================== 
SUBROUTINE initMacro(rho,u,uSqr,image) 
    USE simParam, ONLY: xdim, ydim, umax 
    USE cellConst, ONLY: incylinder,wall 
    implicit none 
    REAL(8), INTENT(INOUT):: rho(ydim,xdim), u(ydim,xdim,0:1), uSqr(ydim,xdim) 
    !REAL(8):: uProf 
    !REAL(8):: randomnumber 
    INTEGER, INTENT(IN):: image(ydim,xdim) 
    INTEGER:: x,y 
    DO x=1,xdim 
    do y=1,ydim 
        IF (image(y,x)/=incylinder .and. image(y,x)/=wall) THEN 
            If (MOD(x, 10)==0 .and. MOD(y, 10)==0) THEN 
            !CALL RANDOM_SEED() 
            !CALL RANDOM_NUMBER(randomnumber) 
                u(y,x,0)=umax*0.5 
            ELSE 
                u(y,x,0)=umax 
            END IF 
            !u(y,x,0)=uProf(y) 
            !u(y,x,0)=umax 
            u(y,x,1)=0.0d0 
        END IF 
    enddo 
    END DO 
    rho =1.0d0 
    uSqr=u(:,:,0)*u(:,:,0)+u(:,:,1)*u(:,:,1) 
END SUBROUTINE initMacro 
 
 
!======================================================== 
! Compute equilibrium distribution 
!======================================================== 
SUBROUTINE computeFeq(fEq,rho,u,uSqr) 
    USE D2Q9COnst, ONLY: t, v 
    USE simParam, ONLY: xdim, ydim 
    implicit none 
    REAL(8), INTENT(IN):: rho(ydim,xdim), uSqr(ydim,xdim), u(ydim,xdim,0:1) 
    REAL(8), INTENT(INOUT):: fEq(ydim,xdim,0:8) 
    INTEGER:: i, x, y 
    REAL(8):: uxy 
    do i = 0, 8 
      do x = 1, xdim 
        do y = 1, ydim 
          uxy=u(y,x,0)*v(i,0)+u(y,x,1)*v(i,1) 
          fEq(y,x,i)=t(i)*rho(y,x)*(1.0d0+3.0d0*uxy+4.5d0*uxy*uxy-1.5d0*uSqr(y,x)) 
        enddo 
      enddo 
    enddo 
END SUBROUTINE computeFeq 
 
 
!======================================================== 
! Compute density and velocity from distribution functions 
!======================================================== 
SUBROUTINE computeMacros(f,rho,u,uSqr) 
    USE simParam, ONLY: xdim, ydim 
 
    implicit none 
    REAL(8), INTENT(IN):: f(ydim,xdim,0:8) 
    REAL(8), INTENT(INOUT):: u(ydim,xdim,0:1), rho(ydim, xdim), uSqr(ydim, xdim) 
 
    INTEGER:: x,y 
    do x = 1, xdim 
      do y = 1, ydim 
        rho(y,x)=f(y,x,0)+f(y,x,1)+f(y,x,2)+f(y,x,3)+f(y,x,4)+f(y,x,5)+f(y,x,6)+f(y,x,7)+f(y,x,8) 
        u(y,x,0)=(f(y,x,1)-f(y,x,3)+f(y,x,5)-f(y,x,6)-f(y,x,7)+f(y,x,8))/rho(y,x) 
        u(y,x,1)=(f(y,x,2)-f(y,x,4)+f(y,x,5)+f(y,x,6)-f(y,x,7)-f(y,x,8))/rho(y,x) 
        uSqr(y,x)=u(y,x,0) * u(y,x,0) + u(y,x,1) * u(y,x,1) 
      enddo 
    enddo 
END SUBROUTINE computeMacros 
 
 
!======================================================== 
! Implement Bounce-back on upper/lower boundaries 
!======================================================== 
SUBROUTINE boundaries(f,image) 
    USE D2Q9Const, ONLY: opposite 
    USE cellConst, ONLY: wall 
    USE simParam, ONLY: xdim, ydim 
    implicit none 
 
    INTEGER, INTENT(IN):: image(ydim,xdim) 
    REAL(8), INTENT(INOUT):: f(ydim,xdim,0:8) 
    REAL(8):: fTmp(0:8) 
    INTEGER:: i, x, y 
 
    do x = 1, xdim 
        do y = 1, ydim 
            if (image(y,x) == wall) then 
                do i = 0, 8 
                    fTmp(i) = f(y,x,opposite(i)) 
                end do 
                do i = 0, 8 
                    f(y,x,i) = fTmp(i) 
                end do 
            end if 
        end do 
    end do 
END SUBROUTINE boundaries 
 
 
!======================================================== 
! Use Zou/He boundary condition to implement Dirichlet 
!       boundaries on inlet/outlet 
!======================================================== 
SUBROUTINE inletOutlet(f,rho,u,image) 
    USE cellConst, ONLY: inlet, outlet 
    USE simParam 
    implicit none 
    REAL(8), INTENT(INOUT):: f(ydim,xdim,0:8),u(ydim,xdim,0:1),rho(ydim,xdim) 
    INTEGER, INTENT(IN):: image(ydim,xdim) 
    REAL(8):: uProf 
    INTEGER:: x,y 
    do x=1,xdim 
      do y=1,ydim 
        if(image(y,x) == inlet) then 
          u(y,x,0) = uProf(y) 
          u(y,x,1) = 0.0d0 
          CALL inletZou(f(y,x,:),u(y,x,:),rho(y,x)) 
        elseif(image(y,x) == outlet) then 
          u(y,x,0) = uProf(y) 
          u(y,x,1) = 0.0d0 
          CALL outletZou(f(y,x,:),u(y,x,:),rho(y,x)) 
        endif 
      enddo 
    enddo 
 
CONTAINS 
 
!======================================================== 
! Zou/He boundary on inlet 
!======================================================== 
    SUBROUTINE inletZou(f,u,rho) 
        implicit none 
        REAL(8), INTENT(INOUT):: f(0:8),rho 
        REAL(8), INTENT(IN):: u(0:1) 
        REAL(8):: fInt, fInt2 
        fInt   = f(0) + f(2) + f(4) 
        fInt2  = f(3) + f(6) + f(7) 
        rho    = (fInt + 2.0d0 * fInt2) / (1.0d0 - u(0)) 
        CALL zouWestWall(f,rho,u) 
    END SUBROUTINE inletZou 
 
    SUBROUTINE zouWestWall(f,rho,u) 
        implicit none 
        REAL(8), INTENT(INOUT):: f(0:8) 
        REAL(8), INTENT(IN):: rho, u(0:1) 
        REAL(8):: fDiff, rhoUx, rhoUy 
        fDiff = 0.5d0 * (f(2) - f(4)) 
        rhoUx = rho * u(0) / 6.0d0 
        rhoUy = 0.5d0 * rho * u(1) 
        f(1) = f(3) + 4.0d0 * rhoUx 
        f(5) = f(7) - fDiff + rhoUx + rhoUy 
        f(8) = f(6) + fDiff + rhoUx - rhoUy 
    END SUBROUTINE zouWestWall 
 
 
!======================================================== 
! Zou/He boundary on outlet 
!======================================================== 
    SUBROUTINE outletZou(f,u,rho) 
        implicit none 
        REAL(8), INTENT(INOUT):: f(0:8),rho,u(0:1) 
        REAL(8):: fInt, fInt2 
        fInt  = f(0) + f(2) + f(4) 
        fInt2 = f(1) + f(8) + f(5) 
        rho   = (fInt + 2.0d0 * fInt2) / (1.0d0 + u(0)) 
        CALL zouEastWall(f,rho,u) 
    END SUBROUTINE outletZou 
 
    SUBROUTINE zouEastWall(f,rho,u) 
        implicit none 
        REAL(8), INTENT(INOUT):: f(0:8) 
        REAL(8), INTENT(IN):: rho, u(0:1) 
        REAL(8):: fDiff, rhoUx, rhoUy 
        fDiff = 0.5d0 * (f(2) - f(4)) 
        rhoUx = rho * u(0) / 6.0d0 
        rhoUy = 0.5d0 * rho * u(1) 
        f(3) = f(1) - 4.0d0 * rhoUx 
        f(7) = f(5) + fDiff - rhoUx - rhoUy 
        f(6) = f(8) - fDiff - rhoUx + rhoUy 
    END SUBROUTINE zouEastWall 
 
END SUBROUTINE inletOutlet 
 
 
!======================================================== 
! Computation of Poiseuille profile for the inlet/outlet 
!======================================================== 
FUNCTION uProf(y) 
    USE simParam, ONLY: ydim, umax 
    implicit none 
    INTEGER, INTENT(IN):: y 
    REAL(8):: radius, uProf 
    radius = dble(ydim-1) * 0.5d0 
    uProf  = -umax * ((abs(1 - dble(y-1) / radius))**2 - 1.0d0) 
END FUNCTION uProf 
 
 
!======================================================== 
! Streaming step: the population functions are shifted 
!       one site along their corresponding lattice direction 
!       (no temporary memory is needed) 
!======================================================== 
SUBROUTINE streaming(f,v,image) 
    USE simParam 
    USE cellConst, ONLY: incylinder 
    implicit none 
    INTEGER::y,x,yd,xd,k 
    REAL(8), INTENT(INOUT):: f(ydim,xdim,0:8) 
    INTEGER, INTENT(IN):: v(0:8,0:1) 
    REAL(8):: f_temp(ydim,xdim,0:8) 
    INTEGER, INTENT(IN):: image(ydim,xdim) 
!	 ------------------------------------- 
    f_temp = f 
    DO k = 0,8 
        DO x = 1,xdim 
                DO y = 1,ydim 
                yd = y - v(k,1) 
                xd = x - v(k,0) 
                IF (yd>=1 .and. yd<=ydim .and. xd>=1 .and. xd<=xdim .and. image(y,x) /= incylinder) THEN 
                    f(y,x,k)=f_temp(yd,xd,k) 
                END IF 
            END DO 
        END DO 
    END DO 
 
END SUBROUTINE streaming 
 
 
!======================================================== 
! LBGK collision step 
!======================================================== 
SUBROUTINE collide(f,fEq,omega,image) 
    USE simParam, ONLY: xdim, ydim 
    USE cellConst, ONLY: incylinder 
    implicit none 
    INTEGER, INTENT(IN):: image(ydim,xdim) 
    REAL(8), INTENT(IN):: fEq(ydim,xdim,0:8), omega 
    REAL(8), INTENT(INOUT):: f(ydim,xdim,0:8) 
    INTEGER:: x,y,i 
    do i=0,8 
      do x=1,xdim 
        do y=1,ydim 
          if(image(y,x) /= incylinder) f(y,x,i)=(1.0d0-omega)*f(y,x,i)+omega*feq(y,x,i) 
        enddo 
      enddo 
    enddo 
END SUBROUTINE collide 
 
 
!======================================================== 
! Write the components of the velocity to a text file, 
!       with indices (x,y) 
!======================================================== 
SUBROUTINE writeOutput(u,tStep,vorticity) 
    USE simParam, ONLY: xdim, ydim, Re1, umax 
    implicit none 
    INTEGER, INTENT(IN):: tStep 
    REAL(8), INTENT(IN):: u(ydim,xdim,0:1) 
    REAL(8), INTENT(IN):: vorticity(ydim,xdim) 
    INTEGER:: x,y 
    character (LEN=100):: fileName 
101 format(2i10,3f20.10) 
    write(fileName, '(I7,A8,F7.1,A8,F7.3,A8,I7,A8,I7)') tStep, " RE = ", Re1,& 
          " umax = ", umax, " xdim = ", xdim, " ydim = ", ydim 
    fileName = adjustl(fileName) 
    open(14,file='outputUx_'//trim(fileName)//'.plt') 
    write(14,*)'TITLE=LBfor' 
    write(14,*)'VARIABLES="x","y","u","v","w"' 
    write(14,*)'ZONE F=POINT, I=',ydim,', J=',xdim 
    do x=1,xdim 
      do y=1,ydim 
        write(14,101) x,y,u(y,x,0),u(y,x,1),vorticity(y,x) 
      enddo 
    enddo 
    close(14) 
END SUBROUTINE writeOutput 
 
 
!======================================================== 
! Write the flow geometry to a file 
!======================================================== 
SUBROUTINE writeImage(image) 
    USE simParam, ONLY: xdim, ydim 
    implicit none 
 
    INTEGER, INTENT(IN):: image(ydim,xdim) 
 
    INTEGER:: x,y 
 
    open(13,file='outputImage.dat') 
    do x=1, xdim 
        do y=1, ydim 
            write(13,102) image(y,x) 
        end do 
    end do 
102     format(3i10) 
    close(15) 
END SUBROUTINE writeImage 
 
 
!======================================================== 
! Print out simulation PARAMETERs to screen 
!======================================================== 
SUBROUTINE writeInput(omega) 
    USE simParam 
    implicit none 
    REAL(8), INTENT(IN):: omega 
    write(*,*) 'xdim                 = ', xdim 
    write(*,*) 'ydim                 = ', ydim 
    write(*,*) 'Obstacle X           = ', obstX 
    write(*,*) 'Obstacle Y           = ', obstY 
    write(*,*) 'Obstacle Radius      = ', obstR 
    write(*,*) 'tmax                 = ', tmax 
    write(*,*) 'umax                 = ', umax 
    write(*,*) 'Re                   = ', Re1 
    write(*,*) 'omega                = ', omega 
END SUBROUTINE writeInput 
 
SUBROUTINE computestreamline(u, streamline,image) 
    USE simParam, ONLY: xdim, ydim, obstY 
    USE cellConst, ONLY: incylinder 
    implicit none 
    REAL(8), INTENT(IN):: u(ydim,xdim,0:1) 
    REAL(8), INTENT(INOUT):: streamline(ydim,xdim) 
    INTEGER, INTENT(IN):: image(ydim,xdim) 
    INTEGER:: x,y 
    streamline(1,1)=0.0d0 
    DO x = 2, xdim 
        streamline(1,x) = streamline(1, x-1) - u(1, x-1, 1) * 1.0 
    END DO 
    DO x = 1, xdim 
        DO y = 2, obstY 
            If(image(y-1,x) /= incylinder) THEN 
            streamline(y,x) = streamline(y-1,x) + u(y-1, x, 0) * 1.0 
            END IF 
        END DO 
    END DO 
    streamline(ydim,1)=0.0d0 
    DO x = 2, xdim 
        streamline(ydim,x) = streamline(ydim, x-1) - u(ydim, x-1, 1) * 1.0 
    END DO 
    DO x = 1, xdim 
        DO y = ydim-1, obstY,-1 
            If(image(y,x) /= incylinder) THEN 
            streamline(y,x) = streamline(y+1,x) + u(y+1, x, 0) * 1.0 
            END IF 
        END DO 
    END DO 
 
 
END SUBROUTINE computestreamline 
 
SUBROUTINE computemeq(meq,rho,u,uSqr) 
    USE simParam, ONLY: xdim, ydim 
    implicit none 
    REAL(8), INTENT(IN):: rho(ydim,xdim), uSqr(ydim,xdim), u(ydim,xdim,0:1) 
    REAL(8), INTENT(INOUT):: meq(ydim,xdim,0:8) 
    INTEGER:: i, x, y 
    do i = 0, 8 
      do x = 1, xdim 
        do y = 1, ydim 
            SELECT CASE(i) 
            CASE(0) 
                meq(y,x,i)=rho(y,x) 
            CASE(1) 
                meq(y,x,i)=rho(y,x)*(-2.0d0+3.0d0*usqr(y,x)) 
            CASE(2) 
                meq(y,x,i)=rho(y,x)*(1.0d0-3.0d0*usqr(y,x)) 
            CASE(3) 
                meq(y,x,i)=rho(y,x)*u(y,x,0) 
            CASE(4) 
                meq(y,x,i)=-rho(y,x)*u(y,x,0) 
            CASE(5) 
                meq(y,x,i)=rho(y,x)*u(y,x,1) 
            CASE(6) 
                meq(y,x,i)=-rho(y,x)*u(y,x,1) 
            CASE(7) 
                meq(y,x,i)=rho(y,x)*(u(y,x,0)*u(y,x,0)-u(y,x,1)*u(y,x,1)) 
            CASE(8) 
                meq(y,x,i)=rho(y,x)*u(y,x,0)*u(y,x,1) 
            CASE DEFAULT 
                meq(y,x,i)=0.0d0 
            END SELECT 
         END DO 
      enddo 
    enddo 
END SUBROUTINE computemeq 
 
SUBROUTINE coll_mrt(f,meq,omega,image) 
    USE simParam, ONLY: xdim, ydim 
    USE cellConst, ONLY: incylinder 
    implicit none 
    INTEGER:: x,y,i 
    INTEGER, INTENT(IN):: image(ydim,xdim) 
    REAL(8), INTENT(IN):: meq(ydim,xdim,0:8), omega 
    REAL(8), INTENT(INOUT):: f(ydim,xdim,0:8) 
    REAL(8):: m(0:8),s(0:8) 
    REAL(8):: d(0:8)=(/9.0d0,36.0d0,36.0d0,6.0d0,12.0d0,6.0d0,12.0d0,4.0d0,4.0d0/) 
 
    s(7:8)=omega 
    s(0)=0.0d0 
    s(1)=1.6d0 
    s(2)=1.8d0 
    s(3:5:2)=0.0d0 
    s(4:6:2)=8*(2-s(7))/(8-s(7)) 
 
      do x=1,xdim 
        do y=1,ydim 
          If(image(y,x) /= incylinder) THEN 
            m(0)=f(y,x,0)+f(y,x,1)+f(y,x,2)+f(y,x,3)+f(y,x,4)& 
                  +f(y,x,5)+f(y,x,6)+f(y,x,7)+f(y,x,8) 
            m(1)=-4.0d0*f(y,x,0)-f(y,x,1)-f(y,x,2)-f(y,x,3)-f(y,x,4)& 
                  +2.0d0*(f(y,x,5)+f(y,x,6)+f(y,x,7)+f(y,x,8)) 
            m(2)=4.0d0*f(y,x,0)-2.0d0*(f(y,x,1)+f(y,x,2)+f(y,x,3)+f(y,x,4))& 
                  +f(y,x,5)+f(y,x,6)+f(y,x,7)+f(y,x,8) 
            m(3)=f(y,x,1)-f(y,x,3)+f(y,x,5)-f(y,x,6)-f(y,x,7)+f(y,x,8) 
            m(4)=-2.0d0*(f(y,x,1)-f(y,x,3))+f(y,x,5)-f(y,x,6)-f(y,x,7)+f(y,x,8) 
            m(5)=f(y,x,2)-f(y,x,4)+f(y,x,5)+f(y,x,6)-f(y,x,7)-f(y,x,8) 
            m(6)=-2.0d0*(f(y,x,2)-f(y,x,4))+f(y,x,5)+f(y,x,6)-f(y,x,7)-f(y,x,8) 
            m(7)=f(y,x,1)-f(y,x,2)+f(y,x,3)-f(y,x,4) 
            m(8)=f(y,x,5)-f(y,x,6)+f(y,x,7)-f(y,x,8) 
            DO i=0,8 
                m(i)=m(i)-s(i)*(m(i)-meq(y,x,i)) 
                m(i)=m(i)/d(i) 
            END DO 
            f(y,x,0)=m(0)-4*(m(1)-m(2)) 
            f(y,x,1)=m(0)-m(1)-2*(m(2)+m(4))+m(3)+m(7) 
            f(y,x,2)=m(0)-m(1)-2*(m(2)+m(6))+m(5)-m(7) 
            f(y,x,3)=m(0)-m(1)-2*(m(2)-m(4))-m(3)+m(7) 
            f(y,x,4)=m(0)-m(1)-2*(m(2)-m(6))-m(5)-m(7) 
            f(y,x,5)=m(0)+m(1)+m(1)+m(2)+m(3)+m(4)+m(5)+m(6)+m(8) 
            f(y,x,6)=m(0)+m(1)+m(1)+m(2)-m(3)-m(4)+m(5)+m(6)-m(8) 
            f(y,x,7)=m(0)+m(1)+m(1)+m(2)-m(3)-m(4)-m(5)-m(6)+m(8) 
            f(y,x,8)=m(0)+m(1)+m(1)+m(2)+m(3)+m(4)-m(5)-m(6)-m(8) 
          END IF 
        enddo 
      enddo 
 
END SUBROUTINE coll_mrt 
 
SUBROUTINE computevorticity(u, vorticity, image) 
    USE simParam, ONLY: xdim, ydim, obstY, obstX, obstR 
    USE cellConst, ONLY: wall 
    implicit none 
    REAL(8), INTENT(IN):: u(ydim,xdim,0:1) 
    REAL(8), INTENT(INOUT):: vorticity(ydim,xdim) 
    INTEGER, INTENT(IN):: image(ydim,xdim) 
    INTEGER:: x,y 
    DO x = 2, xdim 
        DO y = 2, obstY-obstR-1 
            vorticity(y,x) = (u(y, x, 1) - u(y, x-1, 1)) / 1.0- (u(y, x, 0) - u(y-1, x, 0)) / 1.0 
 
        END DO 
    END DO 
    DO x = 2, xdim 
        DO y = obstY-obstR-1, obstY+obstR+1 
            If(image(y,x) /= wall) THEN 
            vorticity(y,x) = (u(y, x, 1) - u(y, x-1, 1)) / 1.0 - (u(y, x, 0) - u(y-1, x, 0)) / 1.0 
            END IF 
        END DO 
    END DO 
    DO x = 2, xdim 
        DO y = obstY+obstR+1, ydim-1 
            vorticity(y,x) = (u(y, x, 1) - u(y, x-1, 1)) / 1.0 - (u(y, x, 0) - u(y-1, x, 0)) / 1.0 
 
        END DO 
    END DO 
 
 
END SUBROUTINE computevorticity