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


Ignore:
Timestamp:
May 14, 2015 2:52:56 PM (7 years ago)
Author:
cmaggio
Comment:

migrated content from ccs wiki

Legend:

Unmodified
Added
Removed
Modified
  • cypress/Programming/StokesFlow/ex32.f90

    v1 v1  
     1= Ex32.f90 =
     2{{{#!fortran
     3!
     4!  Large Scale Computing
     5!  Stokes Flow in a Cavity
     6!  ex32.f90
     7!
     8!  Need Lapack and Blas Libraries
     9!
     10program StokesFlowInCavity
     11  use stokeslet2d
     12  implicit none
     13  double precision :: t_cost
     14  double precision,parameter :: EPSILON = 0.005 ! blob size
     15  integer,parameter :: INTGRID = 51  ! # of x-grid lines for internal velocity
     16  double precision :: dp,dm
     17  integer :: numpdepth
     18  integer :: numpwidth
     19  integer :: numparticles
     20  double precision,dimension(:),allocatable :: loc ! particle location
     21  double precision,dimension(:),allocatable :: vel ! particle velocity
     22  double precision,dimension(:),allocatable :: foc ! particle force
     23  double precision,dimension(:),allocatable :: mat ! 2Nx2N dense matrix
     24  double precision,dimension(:),allocatable :: cloc ! array for internal points
     25  double precision,dimension(:),allocatable :: cvel; ! array for internal velocity
     26  integer :: i,j
     27  integer nyg
     28  character(len=32) :: arg
     29 
     30  ! Newer version of Fortran can get commadline arguments
     31  call get_command_argument(1,arg)
     32  if (len_trim(arg) == 0) then
     33     call get_command_argument(0,arg)
     34     print *,arg," [Depth of Cavity]"
     35     call exit(-1)
     36  end if
     37 
     38  ! get inputed depth
     39  read(arg,*) dp
     40 
     41  ! # of particles in depth
     42  numpdepth = int(dp / EPSILON + 0.5)
     43 
     44  ! # of particles in width
     45  numpwidth = int(1.0 / EPSILON + 0.5)
     46 
     47  ! total # od particles
     48  numparticles = numpdepth * 2 + numpwidth * 2
     49  print *,"Total # of Particles=",numparticles
     50 
     51  ! Allocate Space
     52  ! DIM=2 defined by 'stokeslet2d.f90'
     53  allocate(loc(numparticles * DIM))
     54  allocate(vel(numparticles * DIM))
     55  allocate(foc(numparticles * DIM))
     56  allocate(mat(numparticles * DIM * numparticles * DIM))
     57 
     58  ! set location & velocity of particles (blob)
     59  do i = 0, numparticles-1
     60    foc(i * DIM + 1) = 0.d0;
     61    foc(i * DIM + 2) = 0.d0;
     62    if ((i >= 0) .AND. (i <= numpwidth)) then ! top
     63       loc(i * DIM + 1) = -0.5 + EPSILON * dble(i) ! x
     64       loc(i * DIM + 2) = 0.d0 !  y
     65       vel(i * DIM + 1) = 1.d0
     66       vel(i * DIM + 2) = 0.d0
     67    else
     68       if (i <= (numpwidth + numpdepth)) then ! right wall
     69          loc(i * DIM + 1) = 0.5 ! x
     70          loc(i * DIM + 2) = -EPSILON * dble(i - numpwidth) ! y
     71          vel(i * DIM + 1) = 0.d0
     72          vel(i * DIM + 2) = 0.d0
     73       else
     74          if (i <= (2 * numpwidth + numpdepth)) then ! bottom
     75             loc(i * DIM + 1) = 0.5 - EPSILON * dble(i - (numpwidth + numpdepth))! x
     76             loc(i * DIM + 2) = -EPSILON * numpdepth ! y
     77             vel(i * DIM + 1) = 0.d0
     78             vel(i * DIM + 2) = 0.d0     
     79          else                 ! left wall
     80             loc(i * DIM + 1) = -0.5 ! x
     81             loc(i * DIM + 2) = -EPSILON * dble((2 * numpwidth + 2 * numpdepth) - i) ! y
     82             vel(i * DIM + 1) = 0.d0
     83             vel(i * DIM + 2) = 0.d0     
     84          end if
     85       end if
     86    end if
     87 end do
     88 
     89 ! make 2Nx2N Matrix
     90 dm=t_cost();
     91 call slet2d_mkMatrix(numparticles,loc,EPSILON,mat)
     92 print *,"setting matrix ",t_cost(),"(sec)"
     93 
     94 ! Sovle linear ststem
     95 call slet2d_solve(numparticles,mat,vel,foc)
     96 !call slet2d_solve_CG(numparticles,mat,vel,foc)
     97 
     98 print *,"Sovle linear ststem",t_cost(),"(sec)"
     99 
     100 ! check the solution
     101 call slet2d_mkMatrix(numparticles,loc,EPSILON,mat)
     102 call check_solution(numparticles,mat,vel,foc)
     103 
     104 ! deallocate big matrix
     105 deallocate(mat)
     106 
     107 ! out particle (blob) data
     108 open(10,FILE='particle.dat')
     109 do i = 0, numparticles-1
     110    write(10,'(6e16.8)') loc(i * DIM + 1),loc(i * DIM + 2), &
     111         vel(i * DIM + 1),vel(i * DIM + 2),  &
     112         foc(i * DIM + 1),foc(i * DIM + 2)
     113 end do
     114 close(10)
     115 
     116 !
     117 !  compute internal velocity
     118 !
     119 nyg = int(EPSILON * dble(numpdepth * (INTGRID - 1)))
     120 !nyg = numpdepth * (INTGRID - 1) + 1
     121 allocate(cvel(INTGRID * nyg * DIM))
     122 allocate(cloc(INTGRID * nyg * DIM))
     123 !print *,"Internal Grid",INTGRID,"x",nyg
     124 
     125 ! setting location
     126 do j = 0, nyg - 1
     127    do i = 0,INTGRID - 1
     128       cloc(DIM * (i + INTGRID * j) + 1) = -0.5 + dble(i) / dble(INTGRID - 1)
     129       cloc(DIM * (i + INTGRID * j) + 2) = -dble(j) / dble(INTGRID - 1)
     130    end do
     131 end do
     132 
     133 ! compute velocities
     134 dm=t_cost();
     135 call slet2d_velocity(numparticles, loc, foc, EPSILON, INTGRID * nyg, cloc, cvel)
     136 print *,"Compute internal velocity",t_cost(),"(sec)"
     137 
     138 ! out velocities
     139 open(20,FILE='res.dat')
     140 do j = 0, nyg - 1
     141    do i = 0, INTGRID - 1
     142       write(20,'(4e16.8)') cloc(DIM * (i + INTGRID * j) + 1), cloc(DIM * (i + INTGRID * j) + 2),&
     143            cvel(DIM * (i + INTGRID * j) + 1), cvel(DIM * (i + INTGRID * j) + 2)
     144    end do
     145 end do
     146 close(20)
     147 
     148 ! free
     149 deallocate(cloc)
     150 deallocate(cvel)
     151 
     152 deallocate(loc)
     153 deallocate(vel)
     154 deallocate(foc)
     155end program StokesFlowInCavity
     156 
     157! for timing
     158double precision function t_cost()
     159  real, save :: ptm
     160  integer,save :: fl = 0
     161  real :: tm,sec;
     162 
     163  call cpu_time(tm);
     164  if (fl == 0) then
     165     ptm = tm
     166  end if
     167  fl = 1;
     168 
     169  sec = tm - ptm
     170 
     171  ptm = tm
     172 
     173  t_cost = sec
     174  return
     175end function t_cost
     176stokeslet2d.f90
     177
     178!
     179!  Large Scale Computing
     180!  Stokeslet for 2D
     181
     182!  Ricardo Cortez,"The Method of Regularized Stokeslets",
     183!  (2001), SIAM J. Sci. Comput., Vol.23, No.4, pp.1204
     184
     185!  assuming mu = 1
     186!
     187module stokeslet2d
     188  implicit none
     189  parameter DIM=2
     190  REAL, PARAMETER :: PI = 3.141592653589793238462643383279502884197169399375
     191contains
     192 
     193  ! make Matrix
     194  subroutine slet2d_mkMatrix(np,loc,ep,mat)
     195    integer :: np
     196    double precision, dimension(:), allocatable :: loc,mat
     197    double precision :: ep
     198    integer :: i,j
     199    double precision :: r
     200    double precision :: dx,dy
     201    double precision tr1,tr2
     202 
     203    !zeros mat
     204    do i = 1,  DIM * DIM * np * np
     205       mat(i) = 0.d0
     206    end do
     207 
     208    ! make mat
     209    ! Because of symmetric, we can loop over a half 
     210    do i = 0, np-1
     211       do j = i, np-1
     212          dx = loc(i * DIM + 1) - loc(j * DIM + 1)
     213          dy = loc(i * DIM + 2) - loc(j * DIM + 2)
     214          r = sqrt(dx * dx + dy * dy);
     215 
     216          tr1 = term1(r,ep) / (4.d0 * PI);
     217          tr2 = term2(r,ep) / (4.d0 * PI);
     218 
     219          ! diagoanl elements and upper half
     220          mat(i * DIM + 1  + j * DIM      * np * DIM) = &
     221               mat(i * DIM + 1  + j * DIM      * np * DIM) - tr1 + tr2 * dx * dx
     222          mat(i * DIM + 1 + (j * DIM + 1) * np * DIM) = &
     223               mat(i * DIM + 1 + (j * DIM + 1) * np * DIM) + tr2 * dx * dy
     224          mat(i * DIM + 2 +  j * DIM      * np * DIM) = &
     225               mat(i * DIM + 2 +  j * DIM      * np * DIM) + tr2 * dx * dy
     226          mat(i * DIM + 2 + (j * DIM + 1) * np * DIM) = &
     227               mat(i * DIM + 2 + (j * DIM + 1) * np * DIM) - tr1 + tr2 * dy * dy
     228          if (i /= j) then
     229             ! lower half
     230             mat(j * DIM + 1 +  i * DIM      * np * DIM) = &
     231                  mat(j * DIM + 1 +  i * DIM      * np * DIM) - tr1 + tr2 * dx * dx
     232             mat(j * DIM + 1 + (i * DIM + 1) * np * DIM) = &
     233                  mat(j * DIM + 1 + (i * DIM + 1) * np * DIM) + tr2 * dx * dy
     234             mat(j * DIM + 2 +  i * DIM      * np * DIM) = &
     235                  mat(j * DIM + 2 +  i * DIM      * np * DIM) + tr2 * dx * dy
     236             mat(j * DIM + 2 + (i * DIM + 1) * np * DIM) = &
     237                  mat(j * DIM + 2 + (i * DIM + 1) * np * DIM) - tr1 + tr2 * dy * dy
     238          end if
     239       end do
     240    end do
     241    return
     242  end subroutine slet2d_mkMatrix
     243 
     244  ! Sovle Liear ststem
     245  ! Use Lapack
     246  subroutine slet2d_solve(np,mat,vel,foc)
     247    external dgesv
     248    integer ::  np
     249    double precision, dimension(:), allocatable :: mat,vel,foc
     250    integer :: nrhs
     251    integer :: n
     252    integer :: lda
     253    integer :: ldb
     254    integer,dimension(:),allocatable :: ipiv;
     255    integer :: i;
     256    integer :: info;
     257 
     258    ! solve Ax=b
     259    ! A (amat) is n x n matrix
     260    nrhs = 1 ! # of RHS is 1
     261    n = DIM * np
     262    lda = DIM * np
     263    ldb = DIM * np
     264    ! The pivot indices that define the permutation matrix P
     265    ! row i of the matrix was interchanged with row IPIV(i).
     266    allocate(ipiv(DIM * np))
     267 
     268    ! set RHS
     269    do i = 1, n
     270       foc(i) = vel(i)
     271    end do
     272 
     273    ! solve with LAPACK
     274    call dgesv(n, nrhs, mat, lda, ipiv, foc, ldb, info)
     275 
     276    ! check
     277    ! if (info == 0) printf("successfully done\n");
     278    if (info < 0) then
     279       print *,"the",-info,"-th argument had an illegal value"
     280       call exit(-1)
     281    end if
     282    if (info > 0) then
     283       print *,"U(",info,",",info,") is exactly zero."
     284       call exit(-1)
     285    end if
     286 
     287    ! free
     288    deallocate(ipiv)
     289    return
     290  end subroutine slet2d_solve
     291 
     292  ! CG use blas
     293  subroutine slet2d_solve_CG(np,mat,b,x)
     294    external dnrm2,dcopy,dgemv,ddot,daxpy
     295    double precision :: ddot, dnrm2
     296    integer :: np
     297    double precision,dimension(:),allocatable :: mat,b,x
     298    double precision,parameter :: tol = 1.d-8
     299    double precision,dimension(:),allocatable :: r
     300    double precision,dimension(:),allocatable :: p
     301    double precision,dimension(:),allocatable :: Ap
     302    double precision :: rr,pAp,rr1,alpha,beta
     303    integer :: i,j,count
     304    double precision :: al,bt,normb
     305    character(len=1) ::  tr = "N"
     306    integer :: inc,n
     307 
     308    n = DIM * np
     309    inc = 1
     310 
     311    ! Check if the solution is trivial.
     312    normb = dnrm2(n, b, inc)
     313    if (normb <= 0.0) then
     314       do i = 1, n
     315          x(i) = 0.d0
     316       end do
     317       return
     318    end if
     319 
     320    ! Allocate Working Spaces
     321    allocate(r(n))
     322    allocate(p(n))
     323    allocate(Ap(n))
     324 
     325    ! compute r0 = b - Ax
     326    al = -1.d0
     327    bt = 1.d0
     328    call dcopy(n,b,inc,r,inc) ! r = b
     329    call dgemv(tr,n,n,al,mat,n,x,inc,bt,r,inc) ! r=r-Ax
     330 
     331    rr = 0.d0
     332    call dcopy(n,r,inc,p,inc) ! p = r
     333    rr = ddot(n,r,inc,r,inc) ! rr = r.r
     334    print *,"rr=",rr
     335 
     336    ! cg iteration
     337    count = 0
     338    do while(rr  > tol * tol * normb * normb)
     339       ! Ap = A*p
     340       al = 1.d0
     341       bt = 0.d0
     342       call dgemv(tr,n,n,al,mat,n,p,inc,bt,Ap,inc)
     343 
     344       ! alpha = r.r / p.Ap
     345       pAp = 0.d0
     346       pAp = ddot(n,p,inc,Ap,inc)  ! pAp = p.Ap
     347       alpha = rr / pAp
     348 
     349       ! Beta
     350       rr1 = 0.d0
     351       call daxpy(n,alpha,p,inc,x,inc) ! x += alpha * p
     352       al = -alpha
     353       call daxpy(n,al,Ap,inc,r,inc) ! r -= alpha * Ap
     354       rr1 = ddot(n,r,inc,r,inc) ! rr1 = r.r
     355 
     356       beta = rr1 / rr
     357       ! p = r + beta * p  no blas routine :(
     358       do i = 1,n
     359          p(i) = r(i) + beta * p(i)
     360       end do
     361 
     362       rr = rr1
     363       count = count + 1
     364    end do
     365 
     366    ! Deallocate Working Spaces
     367    deallocate(r)
     368    deallocate(p)
     369    deallocate(Ap)
     370    return
     371  end subroutine slet2d_solve_CG
     372 
     373  subroutine check_solution(np,mat,b,x)
     374    external dcopy,dgemv,dnrm2
     375    double precision :: dnrm2
     376    integer :: np
     377    double precision,dimension(:),allocatable :: mat,b,x
     378    double precision,dimension(:),allocatable :: r
     379    double precision :: al,bt,rr
     380    character(len=1) :: tr = "N"
     381    integer :: inc,n
     382 
     383    n = DIM * np
     384    inc = 1
     385 
     386    ! Allocate Working Spaces
     387    allocate(r(n))
     388 
     389    ! compute r0 = b - Ax
     390    al = -1.d0
     391    bt = 1.d0
     392    call dcopy(n,b,inc,r,inc) ! r = b
     393    call dgemv(tr,n,n,al,mat,n,x,inc,bt,r,inc)
     394 
     395    ! rr = |r|
     396    rr = dnrm2(n, r, inc)
     397 
     398    print *,"|b - Ax|=",rr
     399    deallocate(r)
     400    return
     401  end subroutine check_solution
     402 
     403  !
     404  ! compute velocity
     405  subroutine slet2d_velocity(np, loc, foc, ep, nump, cloc,cvel)
     406    integer :: np
     407    double precision,dimension(:),allocatable :: loc, foc
     408    double precision :: ep
     409    integer :: nump
     410    double precision,dimension(:),allocatable :: cloc,cvel
     411    integer :: i,p
     412    double precision :: r
     413    double precision :: dx,dy
     414    double precision :: tr1,tr2
     415 
     416    ! loop for grid
     417    do p = 0, nump-1
     418       ! zeros
     419       cvel(p * DIM + 1) = 0.d0
     420       cvel(p * DIM + 2) = 0.d0
     421 
     422       ! loop for partilces (blob)
     423       do i = 0, np-1
     424          dx = cloc(p * DIM + 1) - loc(i * DIM + 1)
     425          dy = cloc(p * DIM + 2) - loc(i * DIM + 2)
     426          r = sqrt(dx * dx + dy * dy)
     427 
     428          tr1 = term1(r,ep) / (4.d0 * PI)
     429          tr2 = term2(r,ep) / (4.d0 * PI)
     430 
     431          tr2 = tr2 * (foc(i * DIM + 1) * dx + foc(i * DIM + 2) * dy)
     432 
     433          cvel(p * DIM + 1) = cvel(p * DIM + 1) - foc(i * DIM + 1) * tr1 + tr2 * dx
     434          cvel(p * DIM + 2) = cvel(p * DIM + 2) - foc(i * DIM + 2) * tr1 + tr2 * dy
     435       end do
     436    end do
     437    return
     438  end subroutine slet2d_velocity
     439 
     440  double precision function term1(r,ep)
     441    double precision :: r,ep
     442    double precision ::sq
     443 
     444    sq = sqrt(r * r + ep * ep)
     445 
     446    term1 = log(sq + ep) - ep * (sq + 2.d0 * ep) / (sq + ep) / sq
     447    return
     448  end function term1
     449 
     450  double precision function term2(r,ep)
     451    double precision :: r,ep
     452    double precision :: sq
     453 
     454    sq = sqrt(r * r + ep * ep)
     455 
     456    term2 = (sq + 2.d0 * ep) / (sq + ep) / (sq + ep) / sq
     457    return
     458  end function term2
     459  !
     460end module stokeslet2d
     461}}}