| 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 | ! |
| 10 | program 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) |
| 155 | end program StokesFlowInCavity |
| 156 | |
| 157 | ! for timing |
| 158 | double 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 |
| 175 | end function t_cost |
| 176 | stokeslet2d.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 | ! |
| 187 | module stokeslet2d |
| 188 | implicit none |
| 189 | parameter DIM=2 |
| 190 | REAL, PARAMETER :: PI = 3.141592653589793238462643383279502884197169399375 |
| 191 | contains |
| 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 | ! |
| 460 | end module stokeslet2d |
| 461 | }}} |