= Ex32.f90 =
{{{#!fortran
!
!  Large Scale Computing
!  Stokes Flow in a Cavity
!  ex32.f90
!
!  Need Lapack and Blas Libraries
!
program StokesFlowInCavity
  use stokeslet2d
  implicit none
  double precision :: t_cost
  double precision,parameter :: EPSILON = 0.005 ! blob size
  integer,parameter :: INTGRID = 51  ! # of x-grid lines for internal velocity 
  double precision :: dp,dm
  integer :: numpdepth
  integer :: numpwidth
  integer :: numparticles
  double precision,dimension(:),allocatable :: loc ! particle location
  double precision,dimension(:),allocatable :: vel ! particle velocity
  double precision,dimension(:),allocatable :: foc ! particle force
  double precision,dimension(:),allocatable :: mat ! 2Nx2N dense matrix
  double precision,dimension(:),allocatable :: cloc ! array for internal points
  double precision,dimension(:),allocatable :: cvel; ! array for internal velocity
  integer :: i,j
  integer nyg
  character(len=32) :: arg
 
  ! Newer version of Fortran can get commadline arguments
  call get_command_argument(1,arg)
  if (len_trim(arg) == 0) then
     call get_command_argument(0,arg)
     print *,arg," [Depth of Cavity]"
     call exit(-1)
  end if
 
  ! get inputed depth
  read(arg,*) dp
 
  ! # of particles in depth 
  numpdepth = int(dp / EPSILON + 0.5)
 
  ! # of particles in width 
  numpwidth = int(1.0 / EPSILON + 0.5)
 
  ! total # od particles
  numparticles = numpdepth * 2 + numpwidth * 2
  print *,"Total # of Particles=",numparticles
 
  ! Allocate Space 
  ! DIM=2 defined by 'stokeslet2d.f90'
  allocate(loc(numparticles * DIM))
  allocate(vel(numparticles * DIM))
  allocate(foc(numparticles * DIM))
  allocate(mat(numparticles * DIM * numparticles * DIM))
 
  ! set location & velocity of particles (blob)
  do i = 0, numparticles-1
    foc(i * DIM + 1) = 0.d0;
    foc(i * DIM + 2) = 0.d0;
    if ((i >= 0) .AND. (i <= numpwidth)) then ! top 
       loc(i * DIM + 1) = -0.5 + EPSILON * dble(i) ! x 
       loc(i * DIM + 2) = 0.d0 !  y 
       vel(i * DIM + 1) = 1.d0
       vel(i * DIM + 2) = 0.d0
    else
       if (i <= (numpwidth + numpdepth)) then ! right wall
          loc(i * DIM + 1) = 0.5 ! x 
          loc(i * DIM + 2) = -EPSILON * dble(i - numpwidth) ! y
          vel(i * DIM + 1) = 0.d0
          vel(i * DIM + 2) = 0.d0
       else
          if (i <= (2 * numpwidth + numpdepth)) then ! bottom
             loc(i * DIM + 1) = 0.5 - EPSILON * dble(i - (numpwidth + numpdepth))! x 
             loc(i * DIM + 2) = -EPSILON * numpdepth ! y 
             vel(i * DIM + 1) = 0.d0
             vel(i * DIM + 2) = 0.d0	  
          else                 ! left wall
             loc(i * DIM + 1) = -0.5 ! x
             loc(i * DIM + 2) = -EPSILON * dble((2 * numpwidth + 2 * numpdepth) - i) ! y
             vel(i * DIM + 1) = 0.d0
             vel(i * DIM + 2) = 0.d0	  
          end if
       end if
    end if
 end do
 
 ! make 2Nx2N Matrix 
 dm=t_cost();
 call slet2d_mkMatrix(numparticles,loc,EPSILON,mat)
 print *,"setting matrix ",t_cost(),"(sec)"
 
 ! Sovle linear ststem 
 call slet2d_solve(numparticles,mat,vel,foc)
 !call slet2d_solve_CG(numparticles,mat,vel,foc)
 
 print *,"Sovle linear ststem",t_cost(),"(sec)"
 
 ! check the solution 
 call slet2d_mkMatrix(numparticles,loc,EPSILON,mat)
 call check_solution(numparticles,mat,vel,foc)
 
 ! deallocate big matrix 
 deallocate(mat)
 
 ! out particle (blob) data 
 open(10,FILE='particle.dat')
 do i = 0, numparticles-1
    write(10,'(6e16.8)') loc(i * DIM + 1),loc(i * DIM + 2), &
         vel(i * DIM + 1),vel(i * DIM + 2),  &
         foc(i * DIM + 1),foc(i * DIM + 2)
 end do
 close(10)
 
 ! 
 !  compute internal velocity 
 !
 nyg = int(EPSILON * dble(numpdepth * (INTGRID - 1)))
 !nyg = numpdepth * (INTGRID - 1) + 1
 allocate(cvel(INTGRID * nyg * DIM))
 allocate(cloc(INTGRID * nyg * DIM))
 !print *,"Internal Grid",INTGRID,"x",nyg
 
 ! setting location 
 do j = 0, nyg - 1
    do i = 0,INTGRID - 1
       cloc(DIM * (i + INTGRID * j) + 1) = -0.5 + dble(i) / dble(INTGRID - 1)
       cloc(DIM * (i + INTGRID * j) + 2) = -dble(j) / dble(INTGRID - 1)
    end do
 end do
 
 ! compute velocities 
 dm=t_cost();
 call slet2d_velocity(numparticles, loc, foc, EPSILON, INTGRID * nyg, cloc, cvel)
 print *,"Compute internal velocity",t_cost(),"(sec)"
 
 ! out velocities 
 open(20,FILE='res.dat')
 do j = 0, nyg - 1
    do i = 0, INTGRID - 1
       write(20,'(4e16.8)') cloc(DIM * (i + INTGRID * j) + 1), cloc(DIM * (i + INTGRID * j) + 2),&
            cvel(DIM * (i + INTGRID * j) + 1), cvel(DIM * (i + INTGRID * j) + 2)
    end do
 end do
 close(20)
 
 ! free 
 deallocate(cloc)
 deallocate(cvel)
 
 deallocate(loc)
 deallocate(vel)
 deallocate(foc)
end program StokesFlowInCavity
 
! for timing 
double precision function t_cost()
  real, save :: ptm
  integer,save :: fl = 0
  real :: tm,sec;
 
  call cpu_time(tm);
  if (fl == 0) then
     ptm = tm
  end if
  fl = 1;
 
  sec = tm - ptm
 
  ptm = tm
 
  t_cost = sec
  return
end function t_cost
stokeslet2d.f90

!
!  Large Scale Computing
!  Stokeslet for 2D
!  
!  Ricardo Cortez,"The Method of Regularized Stokeslets", 
!  (2001), SIAM J. Sci. Comput., Vol.23, No.4, pp.1204
!  
!  assuming mu = 1
!
module stokeslet2d
  implicit none
  parameter DIM=2
  REAL, PARAMETER :: PI = 3.141592653589793238462643383279502884197169399375
contains
 
  ! make Matrix 
  subroutine slet2d_mkMatrix(np,loc,ep,mat)
    integer :: np
    double precision, dimension(:), allocatable :: loc,mat
    double precision :: ep
    integer :: i,j
    double precision :: r
    double precision :: dx,dy
    double precision tr1,tr2
 
    !zeros mat
    do i = 1,  DIM * DIM * np * np
       mat(i) = 0.d0
    end do
 
    ! make mat
    ! Because of symmetric, we can loop over a half  
    do i = 0, np-1
       do j = i, np-1
          dx = loc(i * DIM + 1) - loc(j * DIM + 1)
          dy = loc(i * DIM + 2) - loc(j * DIM + 2)
          r = sqrt(dx * dx + dy * dy);
 
          tr1 = term1(r,ep) / (4.d0 * PI);
          tr2 = term2(r,ep) / (4.d0 * PI); 
 
          ! diagoanl elements and upper half 
          mat(i * DIM + 1  + j * DIM      * np * DIM) = &
               mat(i * DIM + 1  + j * DIM      * np * DIM) - tr1 + tr2 * dx * dx
          mat(i * DIM + 1 + (j * DIM + 1) * np * DIM) = &
               mat(i * DIM + 1 + (j * DIM + 1) * np * DIM) + tr2 * dx * dy
          mat(i * DIM + 2 +  j * DIM      * np * DIM) = &
               mat(i * DIM + 2 +  j * DIM      * np * DIM) + tr2 * dx * dy
          mat(i * DIM + 2 + (j * DIM + 1) * np * DIM) = &
               mat(i * DIM + 2 + (j * DIM + 1) * np * DIM) - tr1 + tr2 * dy * dy
          if (i /= j) then
             ! lower half 
             mat(j * DIM + 1 +  i * DIM      * np * DIM) = &
                  mat(j * DIM + 1 +  i * DIM      * np * DIM) - tr1 + tr2 * dx * dx
             mat(j * DIM + 1 + (i * DIM + 1) * np * DIM) = &
                  mat(j * DIM + 1 + (i * DIM + 1) * np * DIM) + tr2 * dx * dy
             mat(j * DIM + 2 +  i * DIM      * np * DIM) = &
                  mat(j * DIM + 2 +  i * DIM      * np * DIM) + tr2 * dx * dy
             mat(j * DIM + 2 + (i * DIM + 1) * np * DIM) = &
                  mat(j * DIM + 2 + (i * DIM + 1) * np * DIM) - tr1 + tr2 * dy * dy
          end if
       end do
    end do
    return
  end subroutine slet2d_mkMatrix
 
  ! Sovle Liear ststem 
  ! Use Lapack 
  subroutine slet2d_solve(np,mat,vel,foc)
    external dgesv
    integer ::  np
    double precision, dimension(:), allocatable :: mat,vel,foc
    integer :: nrhs
    integer :: n
    integer :: lda
    integer :: ldb
    integer,dimension(:),allocatable :: ipiv;
    integer :: i;
    integer :: info;
 
    ! solve Ax=b 
    ! A (amat) is n x n matrix 
    nrhs = 1 ! # of RHS is 1 
    n = DIM * np
    lda = DIM * np
    ldb = DIM * np
    ! The pivot indices that define the permutation matrix P
    ! row i of the matrix was interchanged with row IPIV(i). 
    allocate(ipiv(DIM * np))
 
    ! set RHS 
    do i = 1, n
       foc(i) = vel(i)
    end do
 
    ! solve with LAPACK 
    call dgesv(n, nrhs, mat, lda, ipiv, foc, ldb, info)
 
    ! check 
    ! if (info == 0) printf("successfully done\n"); 
    if (info < 0) then
       print *,"the",-info,"-th argument had an illegal value"
       call exit(-1)
    end if
    if (info > 0) then
       print *,"U(",info,",",info,") is exactly zero."
       call exit(-1)
    end if
 
    ! free 
    deallocate(ipiv)
    return
  end subroutine slet2d_solve
 
  ! CG use blas 
  subroutine slet2d_solve_CG(np,mat,b,x)
    external dnrm2,dcopy,dgemv,ddot,daxpy
    double precision :: ddot, dnrm2
    integer :: np
    double precision,dimension(:),allocatable :: mat,b,x
    double precision,parameter :: tol = 1.d-8
    double precision,dimension(:),allocatable :: r
    double precision,dimension(:),allocatable :: p
    double precision,dimension(:),allocatable :: Ap
    double precision :: rr,pAp,rr1,alpha,beta
    integer :: i,j,count
    double precision :: al,bt,normb
    character(len=1) ::  tr = "N"
    integer :: inc,n
 
    n = DIM * np
    inc = 1
 
    ! Check if the solution is trivial.
    normb = dnrm2(n, b, inc)
    if (normb <= 0.0) then
       do i = 1, n
          x(i) = 0.d0
       end do
       return
    end if
 
    ! Allocate Working Spaces 
    allocate(r(n))
    allocate(p(n))
    allocate(Ap(n))
 
    ! compute r0 = b - Ax 
    al = -1.d0
    bt = 1.d0
    call dcopy(n,b,inc,r,inc) ! r = b 
    call dgemv(tr,n,n,al,mat,n,x,inc,bt,r,inc) ! r=r-Ax
 
    rr = 0.d0
    call dcopy(n,r,inc,p,inc) ! p = r 
    rr = ddot(n,r,inc,r,inc) ! rr = r.r
    print *,"rr=",rr
 
    ! cg iteration 
    count = 0
    do while(rr  > tol * tol * normb * normb) 
       ! Ap = A*p
       al = 1.d0
       bt = 0.d0
       call dgemv(tr,n,n,al,mat,n,p,inc,bt,Ap,inc)
 
       ! alpha = r.r / p.Ap 
       pAp = 0.d0
       pAp = ddot(n,p,inc,Ap,inc)  ! pAp = p.Ap 
       alpha = rr / pAp
 
       ! Beta 
       rr1 = 0.d0
       call daxpy(n,alpha,p,inc,x,inc) ! x += alpha * p
       al = -alpha
       call daxpy(n,al,Ap,inc,r,inc) ! r -= alpha * Ap
       rr1 = ddot(n,r,inc,r,inc) ! rr1 = r.r
 
       beta = rr1 / rr
       ! p = r + beta * p  no blas routine :( 
       do i = 1,n
          p(i) = r(i) + beta * p(i)
       end do
 
       rr = rr1
       count = count + 1
    end do
 
    ! Deallocate Working Spaces 
    deallocate(r)
    deallocate(p)
    deallocate(Ap)
    return
  end subroutine slet2d_solve_CG
 
  subroutine check_solution(np,mat,b,x)
    external dcopy,dgemv,dnrm2
    double precision :: dnrm2
    integer :: np
    double precision,dimension(:),allocatable :: mat,b,x
    double precision,dimension(:),allocatable :: r
    double precision :: al,bt,rr
    character(len=1) :: tr = "N"
    integer :: inc,n
 
    n = DIM * np
    inc = 1
 
    ! Allocate Working Spaces 
    allocate(r(n))
 
    ! compute r0 = b - Ax 
    al = -1.d0
    bt = 1.d0
    call dcopy(n,b,inc,r,inc) ! r = b 
    call dgemv(tr,n,n,al,mat,n,x,inc,bt,r,inc)
 
    ! rr = |r| 
    rr = dnrm2(n, r, inc)
 
    print *,"|b - Ax|=",rr
    deallocate(r)
    return
  end subroutine check_solution
 
  !
  ! compute velocity 
  subroutine slet2d_velocity(np, loc, foc, ep, nump, cloc,cvel)
    integer :: np
    double precision,dimension(:),allocatable :: loc, foc
    double precision :: ep
    integer :: nump
    double precision,dimension(:),allocatable :: cloc,cvel
    integer :: i,p
    double precision :: r
    double precision :: dx,dy
    double precision :: tr1,tr2
 
    ! loop for grid 
    do p = 0, nump-1
       ! zeros 
       cvel(p * DIM + 1) = 0.d0
       cvel(p * DIM + 2) = 0.d0
 
       ! loop for partilces (blob)
       do i = 0, np-1
          dx = cloc(p * DIM + 1) - loc(i * DIM + 1)
          dy = cloc(p * DIM + 2) - loc(i * DIM + 2)
          r = sqrt(dx * dx + dy * dy)
 
          tr1 = term1(r,ep) / (4.d0 * PI)
          tr2 = term2(r,ep) / (4.d0 * PI)
 
          tr2 = tr2 * (foc(i * DIM + 1) * dx + foc(i * DIM + 2) * dy)
 
          cvel(p * DIM + 1) = cvel(p * DIM + 1) - foc(i * DIM + 1) * tr1 + tr2 * dx
          cvel(p * DIM + 2) = cvel(p * DIM + 2) - foc(i * DIM + 2) * tr1 + tr2 * dy
       end do
    end do
    return
  end subroutine slet2d_velocity
 
  double precision function term1(r,ep)
    double precision :: r,ep
    double precision ::sq
 
    sq = sqrt(r * r + ep * ep)
 
    term1 = log(sq + ep) - ep * (sq + 2.d0 * ep) / (sq + ep) / sq
    return
  end function term1
 
  double precision function term2(r,ep)
    double precision :: r,ep
    double precision :: sq
 
    sq = sqrt(r * r + ep * ep)
 
    term2 = (sq + 2.d0 * ep) / (sq + ep) / (sq + ep) / sq
    return
  end function term2
  !
end module stokeslet2d
}}}