Changes between Initial Version and Version 1 of cypress/Programming/StokesFlow/ex32.f90

Timestamp:

05/14/15 14:52:56 (11 years ago)

Author:

cmaggio

Comment:

migrated content from ccs wiki

cypress/Programming/StokesFlow/ex32.f90

@@ +1 @@
+= 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
+}}}