! Lattice Boltzmann sample, written in Fortran 90
!
! Copyright (C) 2006 Orestis Malaspinas
! Address: EPFL STI ISE LIN, ME A2 398, 1015 Lausanne, Switzerland
! E-mail: orestis.malaspinas@epfl.ch
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License
! as published by the Free Software Foundation; either version 2
! of the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public 
! License along with this program; if not, write to the Free 
! Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
! Boston, MA  02110-1301, USA.

! unsteady.f90:
! This example examines an unsteady flow past a cylinder placed in a channel.
! The cylinder is offset somewhat from the center of the flow to make the
! steady-state symmetrical flow unstable. At the inlet and outlet, a Poiseuille
! profile is imposed on the velocity. At Reynolds numbers around 100,
! an unstable periodic pattern is created, the Karman vortex street.
! Note that with the implemented Zou/He boundary condition, you must
! increase the resolution to keep the simulation stable if you increase
! the Reynolds number.


!	 ========================================================
!	 Constants that identify different cell-types according
!        to the dynamics they implement
!	 ========================================================
MODULE cellConst
    integer, parameter:: fluid = 0, wall = 1, inlet = 10, outlet = 11
END MODULE cellConst


!	 ========================================================
!	 Lattice constants for the D2Q9 lattice
!	 ========================================================
MODULE D2Q9Const
!	 D2Q9 Weights
    double precision,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 = 250
    integer, parameter:: yDim = 50
    integer, parameter:: obstX = xDim/5
    integer, parameter:: obstY = yDim/2
    integer, parameter:: obstR = yDim/10+1

    integer, parameter:: tMax = 1000

    double precision, parameter:: uMax = 0.02d0
    double precision, parameter:: Re = 10.0d0
END MODULE simParam


!	 ========================================================
!	 The main program, implementing a flow past a cylinder
!	 ========================================================

PROGRAM unsteady
    USE simParam, ONLY: xDim, yDim, tMax
    implicit none

    double precision:: omega, time1, time2, timeTot
    double precision, dimension(:,:,:), allocatable:: f, fEq, u
    double precision, dimension(:,:), allocatable:: rho, uSqr
    integer, dimension(:,:), allocatable:: image
    integer:: tStep

    allocate(f(yDim,xDim,0:8))
    allocate(fEq(yDim,xDim,0:8))
    allocate(u(yDim,xDim,0:1))
    allocate(uSqr(yDim,xDim))
    allocate(rho(yDim,xDim))
    allocate(image(yDim,xDim))

    CALL constructImage(image)
    CALL computeOmega(omega)
    CALL writeInput(omega)
    CALL initMacro(rho,u,uSqr)
    CALL computeFeq(fEq,rho,u,uSqr)

    f = fEq

    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 stream(f)
        CALL CPU_TIME(time2)
        timeTot = timeTot + (time2-time1)
    end do

    CALL writeImage(image)
    CALL writeOutput(u,0)
    write(*,*) dble(tMax) * (dble(yDim * xDim)) / timeTot ,'cells per second'

    deallocate(f)
    deallocate(fEq)
    deallocate(u)
    deallocate(uSqr)
    deallocate(rho)
    deallocate(image)
END PROGRAM unsteady


!	 ========================================================
!	 Compute the relaxation parameter from the Reynolds number
!	 ========================================================
SUBROUTINE computeOmega(omega)
    USE simParam, ONLY: Re,uMax,obstR

    implicit none

    double precision, INTENT(INOUT):: omega
    double precision:: nu

    nu    =  uMax * 2.0d0 * dble(obstR) / Re
    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) = wall
        end do
    end do

END SUBROUTINE constructImage


!	 ========================================================
!	 Initialize the simulation to Poiseuille profile at
!        an equilibrium distribution
!	 ========================================================
SUBROUTINE initMacro(rho,u,uSqr)
    USE simParam, ONLY: xDim, yDim

    implicit none

    double precision, INTENT(INOUT):: rho(yDim,xDim), u(yDim,xDim,0:1), uSqr(yDim,xDim)
    double precision:: uProf
    integer:: y

    do y = 1, yDim
        u(y,:,0) = uProf(y)
        u(y,:,1) = 0.0d0
    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

    double precision, INTENT(IN):: rho(yDim,xDim), uSqr(yDim,xDim), u(yDim,xDim,0:1)
    double precision, INTENT(INOUT):: fEq(yDim,xDim,0:8)
    integer:: i, x, y
    double precision:: 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))
            end do
        end do
    end do
END SUBROUTINE computeFeq


!	 ========================================================
!	 Compute density and velocity from distribution functions
!	 ========================================================
SUBROUTINE computeMacros(f,rho,u,uSqr)
    USE simParam, ONLY: xDIm, yDim
    implicit none

    double precision, INTENT(IN):: f(yDim,xDim,0:8)
    double precision, 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)
        end do
    end do
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)
    double precision, INTENT(INOUT):: f(yDim,xDim,0:8)
    double precision:: 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

    double precision, INTENT(INOUT):: f(yDim,xDim,0:8), u(yDim,xDim,0:1), rho(yDim,xDim)
    integer, INTENT(IN):: image(yDim,xDim)

    double precision:: 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))
            else if (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))
            end if
        end do
    end do

CONTAINS


    !	 ========================================================
    !	 Zou/He boundary on inlet
    !	 ========================================================
    SUBROUTINE inletZou(f,u,rho)
        implicit none

        double precision, INTENT(INOUT):: f(0:8),rho
        double precision, INTENT(IN):: u(0:1)
        double precision:: 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

        double precision, INTENT(INOUT):: f(0:8)
        double precision, INTENT(IN):: rho, u(0:1)
        double precision:: 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

        double precision, INTENT(INOUT):: f(0:8),rho,u(0:1)
        double precision:: 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

        double precision, INTENT(INOUT):: f(0:8)
        double precision, INTENT(IN):: rho, u(0:1)
        double precision:: 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
    double precision:: 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 stream(f)
    USE simParam
    implicit none

    double precision, INTENT(INOUT):: f(yDim,xDim,0:8)
    double precision:: periodicHor(yDim), periodicVert(xDim)

!	 -------------------------------------
!	 right direction
    periodicHor   = f(:,xDim,1)
    f(:,2:xDim,1) = f(:,1:xDim-1,1)
    f(:,1,1)      = periodicHor
!	 -------------------------------------
!	 up direction
    periodicVert    = f(1,:,2)
    f(1:yDim-1,:,2) = f(2:yDim,:,2)
    f(yDim,:,2)     = periodicVert
!	 -------------------------------------
!	 left direction
    periodicHor     = f(:,1,3)
    f(:,1:xDim-1,3) = f(:,2:xDim,3)
    f(:,xDim,3)     = periodicHor
!	 -------------------------------------
!	 down direction
    periodicVert  = f(yDim,:,4)
    f(2:yDim,:,4) = f(1:yDim-1,:,4)
    f(1,:,4)      = periodicVert
!	 -------------------------------------
!	 up-right direction
    periodicVert = f(1,:,5)
    periodicHor  = f(:,xDim,5)
    f(1:yDim-1,2:xDim,5) = f(2:yDim,1:xDim-1,5)
    f(yDim,2:xDim,5)     = periodicVert(1:xDim-1)
    f(yDim,1,5)          = periodicVert(xDim)
    f(1:yDim-1,1,5)      = periodicHor(2:yDim)
!	 -------------------------------------
!	 up-left direction
    periodicVert = f(1,:,6)
    periodicHor  = f(:,1,6)
    f(1:yDim-1,1:xDim-1,6) = f(2:yDim,2:xDim,6)
    f(yDim,1:xDim-1,6)     = periodicVert(2:xDim)
    f(yDim,xDim,6)         = periodicVert(1)
    f(1:yDim-1,xDim,6)     = periodicHor(2:yDim)
!	 -------------------------------------
!	 down-left direction
    periodicVert = f(yDim,:,7)
    periodicHor  = f(:,1,7)
    f(2:yDim,1:xDim-1,7) = f(1:yDim-1,2:xDim,7)
    f(1,1:xDim-1,7)      = periodicVert(2:xDim)
    f(1,xDim,7)          = periodicVert(1)
    f(2:yDim,xDim,7)     = periodicHor(1:yDim-1)
!	 -------------------------------------
!	 down-right direction
    periodicVert = f(yDim,:,8)
    periodicHor  = f(:,xDim,8)
    f(2:yDim,2:xDim,8) = f(1:yDim-1,1:xDim-1,8)
    f(1,2:xDim,8)      = periodicVert(1:xDim-1)
    f(1,1,8)           = periodicVert(xDim)
    f(2:yDim,1,8)      = periodicHor(1:yDim-1)
END SUBROUTINE stream


!	 ========================================================
!	 LBGK collision step
!	 ========================================================
SUBROUTINE collide(f,fEq,omega,image)
    USE simParam, ONLY: xDim, yDim
    USE cellConst, ONLY: wall
    implicit none

    integer, INTENT(IN):: image(yDim,xDim)
    double precision, INTENT(IN):: fEq(yDim,xDim,0:8), omega
    double precision, 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) /= wall) f(y,x,i) = (1.0d0 - omega) * f(y,x,i) + omega * feq(y,x,i)
            end do
        end do
    end do
END SUBROUTINE collide


!	 ========================================================
!	 Write the components of the velocity to a text file, 
!        with indices (x,y)
!	 ========================================================
SUBROUTINE writeOutput(u,tStep)
    USE simParam, ONLY: xDim, yDim
    implicit none

    integer, INTENT(IN):: tStep
    double precision, INTENT(IN):: u(yDim,xDim,0:1)

    integer:: x,y
    character (LEN=100):: fileName

    write(fileName,*) tStep

    fileName = adjustl(fileName)

    open(14,file='outputUx_'//trim(fileName)//'.dat')
    open(15,file='outputUy_'//trim(fileName)//'.dat')
    do x=1, xDim
        do y=1, yDim
            write(14,101) x,y,u(y,x,0)
            write(15,101) x,y,u(y,x,1)
        end do
    end do
101     format(2i10,f20.10)
    close(14)
    close(15)
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

    double precision, 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                   = ', Re
    write(*,*) 'omega                = ', omega
END SUBROUTINE writeInput
