mac.f90

Go to the documentation of this file.

module dim

    implicit none
    
    integer, parameter  :: nx = 90
    integer, parameter  :: ny = 90
    
    double precision, parameter     :: dx = 1.0d0/(nx-1.0d0)
    double precision, parameter     :: dy = 1.0d0/(ny-1.0d0)
    
    double precision, parameter     :: wconst = 1.0d0/(4.0d0*dx)
    
    integer             :: t
    double precision    :: dt
    
    double precision    :: re
    double precision    :: r
end module dim

module solver
    use dim
    implicit none
    double precision, dimension(nx,ny)  :: u
    double precision, dimension(nx,ny)  :: v
    double precision, dimension(nx,ny)  :: p
    double precision, dimension(nx,ny)  :: Fn
    double precision, dimension(nx,ny)  :: Gn
    double precision, dimension(nx,ny)  :: Q
    
    contains
    
    subroutine initzero()
        u(:,:)  = 0.0d0
        v(:,:)  = 0.0d0
        p(:,:)  = 10.0d0
        fn(:,:) = 0.0d0
        gn(:,:) = 0.0d0
        q(:,:)  = 0.0d0
    end subroutine initzero
    
    subroutine returnu(array)
        double precision, dimension(nx,ny), intent(out)   :: array
        array(:,:) = u(:,:)
    end subroutine returnu
    
    subroutine returnv(array)
        double precision, dimension(nx,ny), intent(out)   :: array
        array(:,:) = v(:,:)
    end subroutine returnv
    
    subroutine returnp(array)
        double precision, dimension(nx,ny), intent(out)   :: array
        array(:,:) = p(:,:)
    end subroutine returnp
    
    subroutine ghostcondition()
        use omp_lib
        use dim
        integer     :: i,j

        !X Boundary condition
        !$omp parallel do shared(u) private(i) schedule(dynamic)
        do j=1,ny
            u(1,j)      = 0.0d0
            u(nx,j)     = 0.0d0
        end do
        !$omp end parallel do

        !Y Boundary condition
        !$omp parallel do shared(v)
        do i=1,nx
            v(i,1)      = 0.0d0
            v(i,ny)     = 0.0d0
        end do
        !$omp end parallel do
    end subroutine ghostcondition
    
    subroutine lidcondition()
        use omp_lib
        use dim
        integer     :: i, j

        !U velocity condition
        !$omp parallel do private(i) shared(u) schedule(dynamic)
        do i=2,nx-1
            u(i,1)      =       - u(i,2)
            u(i,ny)     = 2.0d0 - u(i,ny-1)
        end do
        !$omp end parallel do

        !V Velocity condition
        !$omp parallel do private(j) shared(v) schedule(dynamic)
        do j=2,ny-1
            v(1,j)      =       - v(2,j)
            v(nx,j)     =       - v(nx-1,j)
        end do
        !$omp end parallel do
    end subroutine lidcondition
    
    subroutine calctstep()
        use omp_lib
        use dim
        implicit none
        double precision                :: umax1(2),umax2,vmax1(2),vmax2,dtc,dtr

        !We use different step sizes depending on how long we've been running
        if(t .gt. 10) then
            !Calculate first stability condition
            dtr = r*dx*dy*re

            !get the maximum value of u and v using worksharing
            !$omp workshare
                umax1 = maxval(u)
                umax2 = maxval(umax1)
                vmax1 = maxval(v)
                vmax2 = maxval(vmax1)
            !$omp end workshare

            !Calculate second stability condition
            dtc = (1.0d0 / r) / (re * (abs(umax2) + abs(vmax2))**2 )

            !use the smaller of the two
            dt = min(dtr,dtc)
        else
            !If first few steps then take a very small step
            dt = (r/25.0d0)*dx*dy*re
        end if
    end subroutine calctstep
    
    subroutine calcfngn()
        use dim
        use omp_lib
        implicit none
        integer         :: i, j, err

        !$omp parallel do private(i,j) schedule(dynamic)
        do i=2,nx-1
            do j=2,ny-1
                fn(i,j) = u(i,j)+dt*(((u(i+1,j)-2*u(i,j)+u(i-1,j))/(re*dx*dx))+((u(i,j-1)               &
                          -2.0d0*u(i,j)+u(i,j+1))/(re*dy*dy))-(((u(i,j)+u(i+1,j))*(u(i,j)+u(i+1,j))         &
                          -(u(i-1,j)+u(i,j))*(u(i-1,j)+u(i,j)))/(4.0d0*dx))-((((u(i,j)+u(i,j+1))*(v(i+1,j)  &
                          +v(i,j)))-((u(i,j-1)+u(i,j))*(v(i+1,j-1)+v(i,j-1))))/(4.0d0*dy)))
                gn(i,j) = v(i,j)+dt*(((v(i+1,j)-2.0d0*v(i,j)+v(i-1,j))/(re*dx*dx))+((v(i,j-1)               &
                          -2.0d0*v(i,j)+v(i,j+1))/(re*dy*dy))-(((v(i,j+1)+v(i,j))*(v(i,j+1)+v(i,j))         &
                          -(v(i,j-1)+v(i,j))*(v(i,j-1)+v(i,j)))/(4.0d0*dy))-((((u(i,j)+u(i,j+1))*(v(i+1,j)  &
                          +v(i,j)))-((u(i-1,j)+u(i-1,j+1))*(v(i,j)+v(i-1,j))))/(4.0d0*dx)))
            end do
        end do
        !$omp end parallel do
    end subroutine calcfngn
    
    subroutine calcqn()
        use dim
        use omp_lib
        integer                                         :: i, j, err

        !$omp parallel do private(i,j) schedule(dynamic)
        do i=2,nx-1
            do j=2,ny-1
                q(i,j) =((fn(i,j)-fn(i-1,j)+gn(i,j)-gn(i,j-1))/(dt*dx))
            end do
        end do
        !$omp end parallel do
    end subroutine calcqn
    
    subroutine calcvel()
        use dim
        use omp_lib
        integer                                         :: i, j, err

        !$omp parallel do private(i,j) schedule(dynamic)
        do i=2,nx-2
            do j=2,ny-1
                u(i,j) = fn(i,j) - ( p(i+1,j) - p(i,j) ) * (dt/dx)
            end do
        end do
        !$omp end parallel do

        !$omp parallel do private(i,j) schedule(dynamic)
        do i=2,nx-1
            do j=2,ny-2
                v(i,j) = gn(i,j) - ( p(i,j+1) - p(i,j) ) * (dt/dy)
            end do
        end do
        !$omp end parallel do

    end subroutine calcvel
    
    subroutine vort(w)
        use dim
        implicit none
        double precision, dimension(nx,ny), intent(out) :: w
        integer                                         :: i,j
        
        w(:,:) = 0.0d0
        do j=2,ny-1
            do i=2,nx-1
                w(i,j) = wconst * (v(i+1,j) - v(i-1,j) - u(i,j+1) + u(i,j-1))
            end do
        end do
    end subroutine vort   
    
    subroutine poisson()
        use dim
        implicit none
        double precision, parameter                     :: rf = 1.60d0
        double precision, parameter                     :: conv = 0.00010d0
        integer, parameter                              :: itmax = 200
        double precision                                :: change, pold, po, ch
        integer                                         :: iter, im, jm, k
        integer                                         :: i,j

        iter = 0
        change = 1.0d0
        do while((change .ge. conv))
            iter = iter + 1
            change = 0.0d0

            do i=2,nx-1
                do j=2,ny-1
                    pold=p(i,j)
                    p(i,j) = 0.250d0*((p(i-1,j)+p(i,j-1)+p(i+1,j)+p(i,j+1))-(q(i,j)*dx*dx))
                    p(i,j) = pold + rf*(p(i,j)-pold);

                    !Calculates change
                    if(pold .ne. 0) ch = abs((p(i,j)-pold)/pold)

                    !Stores greatest change value
                    if(ch .gt. change) then
                        change = ch
                        im = i
                        jm = j
                        po = pold
                    end if
                end do
            end do

            !update boundaries on poisson solver
            do k=1,nx
                p(k,1)      = p(k,2)       - ((2.0 * v(k,2))       / (re*dx))
                p(k,ny)     = p(k,ny-1) + ((2.0 * v(k,ny-2)) / (re*dx))
                p(1,k)      = p(2,k)       - ((2.0 * u(2,k))       / (re*dx))
                p(nx,k)     = p(nx-1,k) + ((2.0 * u(nx-2,k)) / (re*dx))
            end do

            if(iter .gt. itmax) then
                goto 100
            end if
        end do
100     continue        
    end subroutine poisson
    
    subroutine parpoisson()
        use omp_lib
        use dim
        implicit none

        double precision, parameter                     :: rf = 1.60d0
        double precision, parameter                     :: conv = 0.00010d0
        integer, parameter                              :: itmax = 500
        double precision                                :: pold, po, ch
        integer                                         :: iter, im, jm, k
        integer                                         :: i,j,tstep
        integer                                         :: rowper, colper, numthreads, istart, iend, myid, kstart, kend
        double precision                                :: change, changethread

        iter = 0
        change = 1.0d0

        !For OpenMP
        numthreads = 4
        call omp_set_num_threads(numthreads)

        !calculate sizes of blocks
        rowper = int(ny/numthreads)
        colper = int(nx/numthreads)

        !Fork the threads
        !$omp parallel private(myid, istart, iend, i, j, k,t, iter, pold, changeThread) shared(change)
        myid = omp_get_thread_num()

        istart = myid*rowper + 1
        if(myid == 0) istart = istart + 1
        iend   = (myid*rowper) + rowper
        if(myid == numthreads-1) iend = iend - 1

        !Main iterative loop
        do while(iter .lt. itmax)
            iter = iter + 1

            !$omp single
            change = 0.0d0
            !$omp end single
            !$omp barrier
            !$omp flush

            changethread = 0.0d0
            do j=istart,iend
                do i=2,nx-1
                    pold=p(i,j)
                    p(i,j) = 0.250d0*((p(i-1,j)+p(i,j-1)+p(i+1,j)+p(i,j+1))-(q(i,j)*dx*dx))
                    p(i,j) = pold + rf*(p(i,j)-pold);

                    !Calculates change
                    if(pold .ne. 0) changethread = max(changethread,abs((p(i,j)-pold)/pold))
                end do
            end do

            !calculates change
            !$omp critical
            change = max(change,changethread)
            !$omp end critical
            !$omp barrier

            if(change .lt. conv) exit

            !$omp barrier

            !update boundaries on poisson solver
            !$omp master
            do k=1,nx
                p(k,1)     = p(k,2)       - ((2.0d0 * v(k,2))       / (re*dx))
                p(k,ny) = p(k,ny-1) + ((2.0d0 * v(k,ny-2)) / (re*dx))
                p(1,k)     = p(2,k)       - ((2.0d0 * u(2,k))       / (re*dx))
                p(nx,k) = p(nx-1,k) + ((2.0d0 * u(nx-2,k)) / (re*dx))
            end do
            !$omp end master

        end do
100     continue
        !$omp barrier
        !$omp end parallel
    end subroutine parpoisson
    
end module solver