1c 2c file muh3.f 3c 4c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5c . . 6c . copyright (c) 1999 by UCAR . 7c . . 8c . University Corporation for Atmospheric Research . 9c . . 10c . all rights reserved . 11c . . 12c . . 13c . MUDPACK version 5.0 . 14c . . 15c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16c 17c 18c ... author and specialist 19c 20c John C. Adams (National Center for Atmospheric Research) 21c email: johnad@ucar.edu, phone: 303-497-1213 22 23c ... For MUDPACK 5.0 information, visit the website: 24c (http://www.scd.ucar.edu/css/software/mudpack) 25c 26c ... purpose (see muh3.d for details) 27c 28c (muh3 is a hybrid multigrid/direct method solver) 29c muh3 attempts to produce a second order finite difference 30c approximation to the three dimensional nonseparable elliptic 31c partial differential equation of the form: 32c 33c cxx(x,y,z)*pxx + cyy(x,y,z)*pyy + czz(x,y,z)*pzz + 34c 35c cx(x,y,z)*px + cy(x,y,z)*py + cz(x,y,z)*pz + 36c 37c ce(x,y,z)*p(x,y,z) = r(x,y,z) 38c 39c 40c ... documentation and test files 41c 42c see the documentation file "muh3.d" for a complete discussion 43c of how to use subroutine muh3. file "tmuh3.f" is a test/driver 44c sample program illustrating use of muh3 45c 46c ... required MUDPACK files 47c 48c mudcom.f, mud3ln.f, mud3pn.f 49c 50c 51 subroutine muh3(iparm,fparm,wk,iwk,coef,bndyc,rhs,phi,mopt,ierror) 52 implicit none 53 integer iwk(*),iparm(23),mopt(4),ierror 54 double precision wk(*),phi(*),rhs(*),fparm(8) 55 integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 56 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 57 +kcycle,iprer,ipost,intpol 58 common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 59 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 60 +kcycle,iprer,ipost,intpol 61 double precision xa,xb,yc,yd,ze,zf,tolmax,relmax 62 common/fmud3/xa,xb,yc,yd,ze,zf,tolmax,relmax 63 integer kpbgn,kcbgn,ktxbgn,ktybgn,ktzbgn,nxk,nyk,nzk,ngrid, 64 + klevel,kcur,kps 65 common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), 66 +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps 67 integer ibeta,ialfa,izmat,idmat 68 common/muh3c/ibeta,ialfa,izmat,idmat 69 integer sxy,sxz,syz,kpspace,kcspace,int 70 integer m,lxy,lxz,lyz,isx,jsy,ksz,iisx,jjsy,kksz 71 integer nx,ny,nz,nxny,itx,ity,itz,k,kb,kkb,kk,ip,ic 72 external coef,bndyc 73 data int / 0 / 74 save int 75 ierror = 1 76 intl = iparm(1) ! set and check intl on ALL calls 77 if (intl*(intl-1).ne.0) return 78 if (int.eq.0) then 79 int = 1 80 if (intl.ne.0) return ! very first call is not intl=0 81 end if 82 ierror = 0 83c 84c set input parameters from iparm,fparm internally 85c 86 intl = iparm(1) 87 nxa = iparm(2) 88 nxb = iparm(3) 89 nyc = iparm(4) 90 nyd = iparm(5) 91 nze = iparm(6) 92 nzf = iparm(7) 93c 94c set grid size params 95c 96 ixp = iparm(8) 97 jyq = iparm(9) 98 kzr = iparm(10) 99 iex = iparm(11) 100 jey = iparm(12) 101 kez = iparm(13) 102c 103c set number of subgrids for mg cycling 104c 105 ngrid = max0(iex,jey,kez) 106 nfx = iparm(14) 107 nfy = iparm(15) 108 nfz = iparm(16) 109 110 iguess = iparm(17) 111 maxcy = iparm(18) 112 method = iparm(19) 113 meth2 = iparm(20) 114 nwork = iparm(21) 115c 116c set floating point params 117c 118 xa = fparm(1) 119 xb = fparm(2) 120 yc = fparm(3) 121 yd = fparm(4) 122 ze = fparm(5) 123 zf = fparm(6) 124 tolmax = fparm(7) 125c 126c set multigrid option parameters 127c 128 kcycle = mopt(1) 129 if (kcycle .eq. 0) then 130c 131c use default settings 132c 133 kcycle = 2 134 iprer = 2 135 ipost = 1 136 intpol = 3 137 else 138 iprer = mopt(2) 139 ipost = mopt(3) 140 intpol = mopt(4) 141 end if 142 if (intl .eq. 0) then ! intialization call 143c 144c check input arguments 145c 146 ierror = 2 ! check boundary condition flags 147 if (max0(nxa,nxb,nyc,nyd,nze,nzf).gt.2) return 148 if (min0(nxa,nxb,nyc,nyd,nze,nzf).lt.0) return 149 if (nxa.eq.0.and.nxb.ne.0) return 150 if (nxa.ne.0.and.nxb.eq.0) return 151 if (nyc.eq.0.and.nyd.ne.0) return 152 if (nyc.ne.0.and.nyd.eq.0) return 153 if (nze.eq.0.and.nzf.ne.0) return 154 if (nze.ne.0.and.nzf.eq.0) return 155 ierror = 3 ! check grid sizes 156 if (ixp.lt.2) return 157 if (jyq.lt.2) return 158 if (kzr.lt.2) return 159c periodic size check 160 if (nxa.eq.0 .and. ixp.lt.3) return 161 if (nyc.eq.0 .and. jyq.lt.3) return 162 if (nze.eq.0 .and. kzr.lt.3) return 163 ierror = 4 164 ngrid = max0(iex,jey,kez) 165 if (iex.lt.1) return 166 if (jey.lt.1) return 167 if (kez.lt.1) return 168 ierror = 13 169 if ((iex-1)*(jey-1)*(kez-1) .eq. 0) then 170 if (maxcy.gt.1) return 171 if (iguess.eq.1) return 172 end if 173 if (ngrid.gt.50) return 174 ierror = 5 175 if (nfx.ne.ixp*2**(iex-1)+1) return 176 if (nfy.ne.jyq*2**(jey-1)+1) return 177 if (nfz.ne.kzr*2**(kez-1)+1) return 178 ierror = 6 179 if (iguess*(iguess-1).ne.0) return 180 ierror = 7 181 if (maxcy.lt.1) return 182 ierror = 8 183 if (method.lt.0 .or. method.gt.10) return 184 if (meth2.lt.0 .or. meth2.gt.3) return 185 ierror = 9 186c 187c compute required work space length 188c 189 m = method 190 isx = 0 191 if ((m-1)*(m-4)*(m-5)*(m-7).eq.0) then 192 isx = 3 193 if (nxa.eq.0) then 194 isx = 5 195 end if 196 end if 197 jsy = 0 198 if ((m-2)*(m-4)*(m-6)*(m-7).eq.0) then 199 jsy = 3 200 if (nyc.eq.0) then 201 jsy = 5 202 end if 203 end if 204 ksz = 0 205 if ((m-3)*(m-5)*(m-6)*(m-7).eq.0) then 206 ksz = 3 207 if (nze.eq.0) then 208 ksz = 5 209 end if 210 end if 211c 212c set scales for planar relaxation 213c 214 iisx = 0 215 jjsy = 0 216 kksz = 0 217 lxy = 0 218 lxz = 0 219 lyz = 0 220 nx = nfx 221 ny = nfy 222 nz = nfz 223 if (method.eq.8) then 224 lxy = 1 225 if (meth2.eq.1.or.meth2.eq.3) then 226 iisx = 3 227 if (nxa.eq.0) iisx = 5 228 end if 229 if (meth2.eq.2.or.meth2.eq.3) then 230 jjsy = 3 231 if (nyc.eq.0) jjsy = 5 232 end if 233 end if 234 if (method.eq.9) then 235 lxz = 1 236 if (meth2.eq.1.or.meth2.eq.3) then 237 iisx = 3 238 if (nxa.eq.0) iisx = 5 239 end if 240 if (meth2.eq.2.or.meth2.eq.3) then 241 kksz = 3 242 if (nze.eq.0) kksz = 5 243 end if 244 end if 245 if (method.eq.10) then 246 lyz = 1 247 if (meth2.eq.1.or.meth2.eq.3) then 248 jjsy = 3 249 if (nyc.eq.0) jjsy = 5 250 end if 251 if (meth2.eq.2.or.meth2.eq.3) then 252 kksz = 3 253 if (nze.eq.0) kksz = 5 254 end if 255 end if 256 sxy = 0 257 sxz = 0 258 syz = 0 259c 260c set subgrid sizes 261c 262 do k=1,ngrid 263 nxk(k) = ixp*2**(max0(k+iex-ngrid,1)-1)+1 264 nyk(k) = jyq*2**(max0(k+jey-ngrid,1)-1)+1 265 nzk(k) = kzr*2**(max0(k+kez-ngrid,1)-1)+1 266 end do 267 kps = 1 268 do kb=1,ngrid 269 k = ngrid-kb+1 270 kpspace = 0 271 kcspace = 0 272 if (method.gt.7) then 273c 274c set spacers for planar relaxation 275c 276 do kkb=2,k 277 kk = k-kkb+1 278 nx = nxk(kk) 279 ny = nyk(kk) 280 nz = nzk(kk) 281 if (method.eq.8) nz = nzk(k) 282 if (method.eq.9) ny = nyk(k) 283 if (method.eq.10) nx = nxk(k) 284 kpspace = kpspace + (nx+2)*(ny+2)*(nz+2) 285 kcspace = kcspace + 8*nx*ny*nz 286 end do 287 end if 288 nx = nxk(k) 289 ny = nyk(k) 290 nz = nzk(k) 291c 292c set pointers 293c 294 kpbgn(k) = kps 295 kcbgn(k) = kpbgn(k)+(nx+2)*(ny+2)*(nz+2)+kpspace 296 ktxbgn(k) = kcbgn(k) + 8*nx*ny*nz + kcspace 297 ktybgn(k) = ktxbgn(k) + isx*nx*ny*nz 298 ktzbgn(k) = ktybgn(k) + jsy*nx*ny*nz 299 kps = ktzbgn(k) + ksz*nx*ny*nz 300c 301c sum space in case planar relaxation in xy or xz or yz 302c 303 sxy = sxy + (6+iisx+jjsy)*nx*ny + (nx+2)*(ny+2) 304 sxz = sxz + (6+iisx+kksz)*nx*nz + (nx+2)*(nz+2) 305 syz = syz + (6+jjsy+kksz)*ny*nz + (ny+2)*(nz+2) 306 end do 307c 308c set pointers for direct at coarse grid 309c 310 nx = ixp+1 311 ny = jyq+1 312 nz = kzr+1 313 ibeta = kps+1 314 nxny = nx*ny 315 if (nze .ne. 0) then 316c nonperiodic z b.c. 317 ialfa = ibeta + nxny*nxny*nz 318 izmat = ialfa 319 idmat = ialfa 320 kps = ialfa+nxny*nxny*nz 321 else 322c periodic z b.c. 323 ialfa = ibeta + nxny*nxny*nz 324 izmat = ialfa+nxny*nxny*nz 325 idmat = izmat+nxny*nxny*(nz-2) 326 kps = idmat+nxny*nxny*(nz-2) 327 end if 328c 329c set and check minimal work space 330c 331 nx = nxk(ngrid) 332 ny = nyk(ngrid) 333 nz = nzk(ngrid) 334 iparm(22)=kps+max0((nx+2)*(ny+2)*(nz+2),lxy*sxy,lxz*sxz,lyz*syz) 335 if (iparm(21) .lt. iparm(22)) return 336 ierror = 10 ! check solution region 337 if (xb.le.xa .or. yd.le.yc .or. zf.le.ze) return 338 ierror = 11 339 if (tolmax .lt. 0.d0) return 340 ierror = 12 ! multigrid parameters 341 if (kcycle.lt.0) return 342 if (min0(iprer,ipost).lt.1) return 343 if ((intpol-1)*(intpol-3).ne.0) return 344 if (max0(kcycle,iprer,ipost).gt.2) then 345 ierror = -5 ! inefficient multigrid cycling 346 end if 347 if (ierror .gt. 0) ierror = 0 ! no fatal errors 348c 349c discretize pde at each grid level 350c 351 do kb=1,ngrid 352 k = ngrid-kb+1 353 nx = nxk(k) 354 ny = nyk(k) 355 nz = nzk(k) 356 ip = kpbgn(k) 357 ic = kcbgn(k) 358 itx = ktxbgn(k) 359 ity = ktybgn(k) 360 itz = ktzbgn(k) 361 klevel = k 362 call dismh3(nx,ny,nz,wk(ic),wk(itx), 363 + wk(ity),wk(itz),bndyc,coef,wk,iwk,ierror) 364 if (method.gt.7) then 365c 366c discretize for planar coarsening 367c 368 do kkb=2,k 369 kk = k-kkb+1 370 ip = ip+(nx+2)*(ny+2)*(nz+2) 371 ic = ic + 8*nx*ny*nz 372 nx = nxk(kk) 373 ny = nyk(kk) 374 nz = nzk(kk) 375 if (method.eq.8) nz = nzk(k) 376 if (method.eq.9) ny = nyk(k) 377 if (method.eq.10) nx = nxk(k) 378 call dismh3(nx,ny,nz,wk(ic),wk(itx), 379 + wk(ity),wk(itz),bndyc,coef,wk,iwk,ierror) 380 end do 381 end if 382 end do 383 return 384 end if ! end of intl=0 initialization call block 385 nx = nfx 386 ny = nfy 387 nz = nfz 388 call muh31(nx,ny,nz,rhs,phi,coef,bndyc,wk,iwk) 389c 390c return if singular pde detected 391c 392 if (ierror .eq. 14) return 393 iparm(23) = itero 394 if (ierror .le. 0) then 395 if (tolmax.gt.0.d0) then 396c 397c set final computed maximum relative difference 398c 399 fparm(8) = relmax 400 if (relmax.gt.tolmax .and. ierror.eq.0) then 401c 402c flag convergence failure 403c 404 ierror = -1 405 end if 406 end if 407 end if 408 return 409 end 410 411 subroutine muh31(nx,ny,nz,rhsf,phif,coef,bndyc,wk,iwk) 412 implicit none 413 integer iwk(*) 414 integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 415 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 416 +kcycle,iprer,ipost,intpol 417 double precision xa,xb,yc,yd,ze,zf,tolmax,relmax 418 integer kpbgn,kcbgn,ktxbgn,ktybgn,ktzbgn,nxk,nyk,nzk,ngrid, 419 + klevel,kcur,kps 420 integer nx,ny,nz,ip,ic,ir,irc,ncx,ncy,ncz,ipc,icc 421 integer i,j,k,kb,jk,kk,kkb,ijk,iter 422 double precision phif(nx,ny,nz),rhsf(nx,ny,nz),wk(*),phmax 423 common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 424 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 425 +kcycle,iprer,ipost,intpol 426 common/fmud3/xa,xb,yc,yd,ze,zf,tolmax,relmax 427 common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), 428 +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps 429 integer ibeta,ialfa,izmat,idmat 430 common/muh3c/ibeta,ialfa,izmat,idmat 431 external coef,bndyc 432 nx = nxk(ngrid) 433 ny = nyk(ngrid) 434 nz = nzk(ngrid) 435 ip = kpbgn(ngrid) 436 ic = kcbgn(ngrid) 437 ir = ic+7*nx*ny*nz 438c 439c set phif,rhsf in wk 440c 441 call swk3(nx,ny,nz,phif,rhsf,wk(ip),wk(ir)) 442 if (iguess.eq.0) then 443c 444c no initial guess at finest grid level! 445c 446 do kb=2,ngrid 447 k = ngrid-kb+1 448 nx = nxk(k+1) 449 ny = nyk(k+1) 450 nz = nzk(k+1) 451 ip = kpbgn(k+1) 452 ic = kcbgn(k+1) 453 ir = kcbgn(k+1)+7*nx*ny*nz 454 ncx = nxk(k) 455 ncy = nyk(k) 456 ncz = nzk(k) 457 ipc = kpbgn(k) 458 icc = kcbgn(k) 459 irc = icc+7*ncx*ncy*ncz 460c 461c transfer down to all grid levels 462c 463 call trsfc3(nx,ny,nz,wk(ip),wk(ir),ncx,ncy,ncz, 464 + wk(ipc),wk(irc)) 465 if (method.gt.7) then 466c 467c transfer down for planar coarsening 468c 469 do kkb=1,k-1 470 kk = k-kkb+1 471 ipc = ip + (nx+2)*(ny+2)*(nz+2) 472 icc = ic + 8*nx*ny*nz 473 ncx = nxk(kk) 474 ncy = nyk(kk) 475 ncz = nzk(kk) 476 if (method.eq.8) then 477 ncz = nzk(k+1) 478 else if (method.eq.9) then 479 ncy = nyk(k+1) 480 else 481 ncx = nxk(k+1) 482 end if 483 irc = icc+7*ncx*ncy*ncz 484 call trsfc3(nx,ny,nz,wk(ip),wk(ir),ncx,ncy,ncz, 485 + wk(ipc),wk(irc)) 486 nx = ncx 487 ny = ncy 488 nz = ncz 489 ip = ipc 490 ir = irc 491 ic = icc 492 end do 493 end if 494 end do 495c 496c adjust right hand side at all grid levels in case 497c rhs or specified b.c. in phi or gbdy changed 498c 499 do kb=1,ngrid 500 k = ngrid-kb+1 501 nx = nxk(k) 502 ny = nyk(k) 503 nz = nzk(k) 504 ip = kpbgn(k) 505 ic = kcbgn(k) 506 call adjmh3(nx,ny,nz,wk(ip),wk(ic),bndyc,coef) 507c 508c adjust for planar grid levels if necessary 509c 510 if (method.gt.7) then 511 do kkb=2,k 512 kk = k-kkb+1 513 ip = ip+(nx+2)*(ny+2)*(nz+2) 514 ic = ic + 8*nx*ny*nz 515 nx = nxk(kk) 516 ny = nyk(kk) 517 nz = nzk(kk) 518 if (method.eq.8) then 519 nz = nzk(k) 520 else if (method.eq.9) then 521 ny = nyk(k) 522 else 523 nx = nxk(k) 524 end if 525 call adjmh3(nx,ny,nz,wk(ip),wk(ic),bndyc,coef) 526 end do 527 end if 528 end do 529c 530c execute one full multigrid cycle 531c 532 do k=1,ngrid-1 533 kcur = k 534 call kcymh3(wk,iwk) 535 nx = nxk(k+1) 536 ny = nyk(k+1) 537 nz = nzk(k+1) 538 ip = kpbgn(k+1) 539 ipc = kpbgn(k) 540 ncx = nxk(k) 541 ncy = nyk(k) 542 ncz = nzk(k) 543c 544c lift or prolong approximation from k to k+1 545c 546 call prolon3(ncx,ncy,ncz,wk(ipc),nx,ny,nz,wk(ip), 547 + nxa,nxb,nyc,nyd,nze,nzf,intpol) 548 end do 549 550 else 551c 552c adjust rhs at finest grid level only 553c 554 nx = nxk(ngrid) 555 ny = nyk(ngrid) 556 nz = nzk(ngrid) 557 ip = kpbgn(ngrid) 558 ic = kcbgn(ngrid) 559 call adjmh3(nx,ny,nz,wk(ip),wk(ic),bndyc,coef) 560 end if 561c 562c execute maxcy more multigrid k cycles from finest level 563c 564 kcur = ngrid 565 do iter=1,maxcy 566 itero = iter 567 call kcymh3(wk,iwk) 568 if (tolmax.gt.0.d0) then 569c 570c error control 571c 572 relmax = 0.d0 573 phmax = 0.d0 574 do k=1,nfz 575 kk = k*(nfx+2)*(nfy+2) 576 do j=1,nfy 577 jk = kk+j*(nfx+2) 578 do i=1,nfx 579 ijk = jk+i+1 580 phmax = max(phmax,abs(wk(ijk))) 581 relmax = max(relmax,abs(wk(ijk)-phif(i,j,k))) 582 phif(i,j,k) = wk(ijk) 583 end do 584 end do 585 end do 586c 587c set maximum relative difference and check for convergence 588c 589 if (phmax.gt.0.d0) relmax = relmax/phmax 590 if (relmax.le.tolmax) return 591 end if 592 end do 593c 594c set final iterate in phif 595c 596 do k=1,nfz 597 kk = k*(nfx+2)*(nfy+2) 598 do j=1,nfy 599 jk = kk+j*(nfx+2) 600 do i=1,nfx 601 ijk = jk+i+1 602 phif(i,j,k) = wk(ijk) 603 end do 604 end do 605 end do 606 return 607 end 608 609 subroutine kcymh3(wk,iwk) 610c 611c perform multigrid k-cycle at kcur level 612c kcycle = 1 corresponds to v cycles 613c kcycle = 2 corresponds to w cycles 614c 615 implicit none 616 double precision wk(*) 617 integer iwk(*) 618 integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 619 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 620 +kcycle,iprer,ipost,intpol 621 common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 622 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 623 +kcycle,iprer,ipost,intpol 624 integer kpbgn,kcbgn,ktxbgn,ktybgn,ktzbgn,nxk,nyk,nzk,ngrid, 625 + klevel,kcur,kps 626 common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), 627 +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps 628 integer ibeta,ialfa,izmat,idmat 629 common/muh3c/ibeta,ialfa,izmat,idmat 630 integer nx,ny,nz,ncx,ncy,ncz,nxny 631 integer kount(50),ip,ic,ipc,irc,nrel,l 632 klevel = kcur 633 if (kcur .eq. 1) then 634 nx = nxk(1) 635 ny = nyk(1) 636 nz = nzk(1) 637 nxny = nx*ny 638 ip = kpbgn(1) 639 ic = kcbgn(1) 640c 641c solve at coarsest grid level with direct method and return 642c 643 if (nze .ne. 0) then 644 call dir3(nxny,nx,ny,nz,wk(ip),wk(ic),wk(ibeta),wk(ialfa), 645 + wk(kps),iwk) 646 return 647 else 648 call dir3p(nxny,nx,ny,nz,wk(ip),wk(ic),wk(ibeta),wk(ialfa), 649 + wk(izmat),wk(idmat),wk(kps),iwk) 650 return 651 end if 652 end if 653c 654c pre-relax at current finest grid level > 1 655c 656 do l = 1,iprer 657 call relmh3(wk) 658 end do 659c 660c restrict residual to kcur-1 661c 662 nx = nxk(klevel) 663 ny = nyk(klevel) 664 nz = nzk(klevel) 665 ip = kpbgn(klevel) 666 ic = kcbgn(klevel) 667 ipc = kpbgn(klevel-1) 668 ncx = nxk(klevel-1) 669 ncy = nyk(klevel-1) 670 ncz = nzk(klevel-1) 671 irc = kcbgn(klevel-1)+7*ncx*ncy*ncz 672c 673c use full weighting with residual restriction 674c 675 call resmh3(nx,ny,nz,wk(ip),wk(ic),ncx,ncy,ncz,wk(ipc),wk(irc), 676 + wk(kps)) 677c 678c set counter for grid levels to zero 679c 680 do l = 1,kcur 681 kount(l) = 0 682 end do 683c 684c set new level and continue k-cycling 685c 686 klevel = kcur-1 687 nrel = iprer 688c 689c kcycle control point 690c 691 1 continue 692c 693c post-relax when kcur revisited 694c 695 if (klevel .eq. kcur) go to 2 696c 697c count "hit" at current level 698c 699 kount(klevel) = kount(klevel)+1 700c 701c relax at current level or use direct method if klevel=1 702c 703 if (klevel .gt. 1) then 704 do l = 1,nrel 705 call relmh3(wk) 706 end do 707 else 708 nx = nxk(1) 709 ny = nyk(1) 710 nz = nzk(1) 711 nxny = nx*ny 712 ip = kpbgn(1) 713 ic = kcbgn(1) 714 if (nze .ne. 0) then 715 call dir3(nxny,nx,ny,nz,wk(ip),wk(ic),wk(ibeta),wk(ialfa), 716 + wk(kps),iwk) 717 else 718 call dir3p(nxny,nx,ny,nz,wk(ip),wk(ic),wk(ibeta),wk(ialfa), 719 + wk(izmat),wk(idmat),wk(kps),iwk) 720 end if 721c 722c insure direct method is not called again at coarse level 723c 724 kount(1) = kcycle+1 725 end if 726 if (kount(klevel) .eq. kcycle+1) then 727c 728c kcycle(iprer,ipost) complete at klevel 729c inject correction to finer grid 730c 731 nx = nxk(klevel+1) 732 ny = nyk(klevel+1) 733 nz = nzk(klevel+1) 734 ip = kpbgn(klevel+1) 735 ncx = nxk(klevel) 736 ncy = nyk(klevel) 737 ncz = nzk(klevel) 738 ipc = kpbgn(klevel) 739 call cor3(nx,ny,nz,wk(ip),ncx,ncy,ncz,wk(ipc), 740 + nxa,nxb,nyc,nyd,nze,nzf,intpol,wk(kps)) 741c 742c reset counter to zero at klevel 743c 744 kount(klevel) = 0 745c 746c ascend to next higher level and set to post-relax there 747c 748 klevel = klevel+1 749 nrel = ipost 750 go to 1 751 else 752c 753c kcycle not complete so descend unless at coarsest 754c 755 if (klevel .gt. 1) then 756 nx = nxk(klevel) 757 ny = nyk(klevel) 758 nz = nzk(klevel) 759 ip = kpbgn(klevel) 760 ic = kcbgn(klevel) 761 ncx = nxk(klevel-1) 762 ncy = nyk(klevel-1) 763 ncz = nzk(klevel-1) 764 irc = kcbgn(klevel-1)+7*ncx*ncy*ncz 765 ipc = kpbgn(klevel-1) 766 call resmh3(nx,ny,nz,wk(ip),wk(ic),ncx,ncy,ncz,wk(ipc), 767 + wk(irc),wk(kps)) 768c 769c pre-relax at next coarser level 770c 771 klevel = klevel-1 772 nrel = iprer 773 go to 1 774 else 775c 776c direct method at coarsest level 777c 778 nx = nxk(1) 779 ny = nyk(1) 780 nz = nzk(1) 781 nxny = nx*ny 782 ip = kpbgn(1) 783 ic = kcbgn(1) 784 if (nze .ne. 0) then 785 call dir3(nxny,nx,ny,nz,wk(ip),wk(ic),wk(ibeta),wk(ialfa), 786 + wk(kps),iwk) 787 else 788 call dir3p(nxny,nx,ny,nz,wk(ip),wk(ic),wk(ibeta),wk(ialfa), 789 + wk(izmat),wk(idmat),wk(kps),iwk) 790 end if 791c 792c inject correction to grid level 2 793c 794 ipc = kpbgn(1) 795 ncx = nxk(1) 796 ncy = nyk(1) 797 ncz = nzk(1) 798 ip = kpbgn(2) 799 nx = nxk(2) 800 ny = nyk(2) 801 nz = nzk(2) 802 call cor3(nx,ny,nz,wk(ip),ncx,ncy,ncz,wk(ipc), 803 + nxa,nxb,nyc,nyd,nze,nzf,intpol,wk(kps)) 804c 805c set to post-relax at level 2 806c 807 nrel = ipost 808 klevel = 2 809 go to 1 810 end if 811 end if 812 2 continue 813c 814c post-relax at kcur level 815c 816 do l = 1,ipost 817 call relmh3(wk) 818 end do 819 return 820 end 821 822 subroutine resmh3(nx,ny,nz,phi,cof,ncx,ncy,ncz,phic,rhsc,resf) 823c 824c compute fully weighted residual restriction in rhsc 825c 826 implicit none 827 integer nx,ny,nz,ncx,ncy,ncz,ic,jc,kc,i,j,k 828 integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 829 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 830 +kcycle,iprer,ipost,intpol 831 common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 832 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 833 +kcycle,iprer,ipost,intpol 834 double precision phi(0:nx+1,0:ny+1,0:nz+1) 835 double precision phic(0:ncx+1,0:ncy+1,0:ncz+1) 836 double precision rhsc(ncx,ncy,ncz),resf(nx,ny,nz),cof(nx,ny,nz,8) 837c 838c initialize phic to zero 839c 840 do kc=0,ncz+1 841 do jc=0,ncy+1 842 do ic=0,ncx+1 843 phic(ic,jc,kc) = 0.d0 844 end do 845 end do 846 end do 847c 848c compute fine grid residual 849c 850!$OMP PARALLEL DO PRIVATE(i,j,k), SHARED(resf,cof,phi,nx,ny,nz) 851 do k=1,nz 852 do j=1,ny 853 do i=1,nx 854 resf(i,j,k) = cof(i,j,k,8)-( 855 + cof(i,j,k,1)*phi(i-1,j,k)+ 856 + cof(i,j,k,2)*phi(i+1,j,k)+ 857 + cof(i,j,k,3)*phi(i,j-1,k)+ 858 + cof(i,j,k,4)*phi(i,j+1,k)+ 859 + cof(i,j,k,5)*phi(i,j,k-1)+ 860 + cof(i,j,k,6)*phi(i,j,k+1)+ 861 + cof(i,j,k,7)*phi(i,j,k)) 862 end do 863 end do 864 end do 865c 866c restrict resf to interior coarse mesh in rhsc 867c using fully weighted residual restriction in 3-d 868c 869 call res3(nx,ny,nz,resf,ncx,ncy,ncz,rhsc,nxa,nxb,nyc,nyd,nze,nzf) 870 return 871 end 872 873 subroutine dismh3(nx,ny,nz,cof,tx,ty,tz,bndyc,coef,wk,iwk,ier) 874c 875c discretize the 3-d elliptic pde 876c 877 implicit none 878 integer nx,ny,nz,ier 879 double precision cof(nx,ny,nz,8) 880 double precision tx(nx,ny,nz,*),ty(ny,nx,nz,*),tz(nz,nx,ny,*), 881 +wk(*) 882 integer iwk(*) 883 integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 884 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 885 +kcycle,iprer,ipost,intpol 886 double precision xa,xb,yc,yd,ze,zf,tolmax,relmax 887 integer kpbgn,kcbgn,ktxbgn,ktybgn,ktzbgn,nxk,nyk,nzk,ngrid, 888 + klevel,kcur,kps 889 common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 890 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 891 +kcycle,iprer,ipost,intpol 892 common/fmud3/xa,xb,yc,yd,ze,zf,tolmax,relmax 893 common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), 894 +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps 895 integer ibeta,ialfa,izmat,idmat 896 common/muh3c/ibeta,ialfa,izmat,idmat 897 double precision dlx,dly,dlz,dlx2,dly2,dlz2,dlxx,dlyy,dlzz,cmin, 898 +cemax,alfmax 899 double precision cxx,cyy,czz,cx,cy,cz,ce,alfa,gbdy,x,y,z,c1,c2,c3, 900 +c4,c5,c6 901 integer i,j,k,l,ist,ifn,jst,jfn,kst,kfn,kbdy 902 integer nxny,nxnz,nynz,im1,jm1,km1 903 external bndyc,coef 904c 905c set current grid increments 906c 907 dlx = (xb-xa)/(nx-1) 908 dlx2 = dlx+dlx 909 dlxx = dlx*dlx 910 dly = (yd-yc)/(ny-1) 911 dly2 = dly+dly 912 dlyy = dly*dly 913 dlz = (zf-ze)/(nz-1) 914 dlz2 = dlz+dlz 915 dlzz = dlz*dlz 916 cmin = 1.0 917 cemax = 0.d0 918c 919c set x,y,z subscript limits to bypass specified boundaries 920c when calling coef or bndyc 921c 922 jst = 1 923 jfn = ny 924 ist = 1 925 ifn = nx 926 kst = 1 927 kfn = nz 928 if (nxa.eq.1) ist = 2 929 if (nxb.eq.1) ifn = nx-1 930 if (nyc.eq.1) jst = 2 931 if (nyd.eq.1) jfn = ny-1 932 if (nze.eq.1) kst = 2 933 if (nzf.eq.1) kfn = nz-1 934 do k=kst,kfn 935 z = ze+(k-1)*dlz 936 do j=jst,jfn 937 y = yc+(j-1)*dly 938 do i=ist,ifn 939 x = xa+(i-1)*dlx 940 call coef(x,y,z,cxx,cyy,czz,cx,cy,cz,ce) 941 cmin = min(cmin,cxx,cyy,czz) 942 cemax = max(abs(ce),cemax) 943c 944c check if pde is "hyperbolic" at finest grid level 945c 946 if (klevel.eq.ngrid) then 947 if ((abs(cx)*dlx .gt. abs(cxx+cxx))) ier = -4 948 if ((abs(cy)*dlx .gt. abs(cyy+cyy))) ier = -4 949 if ((abs(cz)*dlx .gt. abs(czz+czz))) ier = -4 950 end if 951c 952c adjust second order coefficients so that pde is not "hyperbolic" 953c this is especially possible on coarser grids if there are non-zero 954c first order terms 955c 956 cxx = max(cxx,abs(cx)*dlx*0.5d0) 957 cyy = max(cyy,abs(cy)*dly*0.5d0) 958 czz = max(czz,abs(cz)*dlz*0.5d0) 959 c1 = cxx/dlxx-cx/dlx2 960 c2 = cxx/dlxx+cx/dlx2 961 c3 = cyy/dlyy-cy/dly2 962 c4 = cyy/dlyy+cy/dly2 963 c5 = czz/dlzz-cz/dlz2 964 c6 = czz/dlzz+cz/dlz2 965 cof(i,j,k,1) = c1 966 cof(i,j,k,2) = c2 967 cof(i,j,k,3) = c3 968 cof(i,j,k,4) = c4 969 cof(i,j,k,5) = c5 970 cof(i,j,k,6) = c6 971 cof(i,j,k,7) = ce-(c1+c2+c3+c4+c5+c6) 972 end do 973 end do 974 end do 975 if (ier .ne. -4) then 976c 977c set nonfatal error flag if ellipticity test fails 978c 979 if (cmin.le.0.d0) ier = -2 980 end if 981 alfmax = 0.d0 982c 983c adjust equation at mixed b.c. 984c 985 if (nxa.eq.2) then 986 kbdy = 1 987 x = xa 988 i = 1 989 do k=kst,kfn 990 z = ze+(k-1)*dlz 991 do j=jst,jfn 992 y = yc+(j-1)*dly 993 call bndyc(kbdy,y,z,alfa,gbdy) 994 alfmax = max(abs(alfa),alfmax) 995 c1 = cof(i,j,k,1) 996 cof(i,j,k,1) = 0.d0 997 cof(i,j,k,2) = cof(i,j,k,2)+c1 998 cof(i,j,k,7) = cof(i,j,k,7)+dlx2*alfa*c1 999 end do 1000 end do 1001 end if 1002 if (nxb.eq.2) then 1003 kbdy = 2 1004 x = xb 1005 i = nx 1006 do k=kst,kfn 1007 z = ze+(k-1)*dlz 1008 do j=jst,jfn 1009 y = yc+(j-1)*dly 1010 call bndyc(kbdy,y,z,alfa,gbdy) 1011 alfmax = max(abs(alfa),alfmax) 1012 c2 = cof(i,j,k,2) 1013 cof(i,j,k,1) = cof(i,j,k,1)+c2 1014 cof(i,j,k,2) = 0.d0 1015 cof(i,j,k,7) = cof(i,j,k,7)-dlx2*alfa*c2 1016 end do 1017 end do 1018 end if 1019 if (nyc.eq.2) then 1020 kbdy = 3 1021 y = yc 1022 j = 1 1023 do k=kst,kfn 1024 z = ze+(k-1)*dlz 1025 do i=ist,ifn 1026 x = xa+(i-1)*dlx 1027 call bndyc(kbdy,x,z,alfa,gbdy) 1028 alfmax = max(abs(alfa),alfmax) 1029 c3 = cof(i,j,k,3) 1030 cof(i,j,k,3) = 0.d0 1031 cof(i,j,k,4) = cof(i,j,k,4)+c3 1032 cof(i,j,k,7) = cof(i,j,k,7)+dly2*alfa*c3 1033 end do 1034 end do 1035 end if 1036 if (nyd.eq.2) then 1037 kbdy = 4 1038 y = yd 1039 j = ny 1040 do k=kst,kfn 1041 z = ze+(k-1)*dlz 1042 do i=ist,ifn 1043 x = xa+(i-1)*dlx 1044 call bndyc(kbdy,x,z,alfa,gbdy) 1045 alfmax = max(abs(alfa),alfmax) 1046 c4 = cof(i,j,k,4) 1047 cof(i,j,k,3) = cof(i,j,k,3)+c4 1048 cof(i,j,k,4) = 0.d0 1049 cof(i,j,k,7) = cof(i,j,k,7)-dly2*c4*alfa 1050 end do 1051 end do 1052 end if 1053 if (nze.eq.2) then 1054 kbdy = 5 1055 z = ze 1056 k = 1 1057 do j=jst,jfn 1058 y = yc+(j-1)*dly 1059 do i=ist,ifn 1060 x = xa+(i-1)*dlx 1061 call bndyc(kbdy,x,y,alfa,gbdy) 1062 alfmax = max(abs(alfa),alfmax) 1063 c5 = cof(i,j,k,5) 1064 cof(i,j,k,5) = 0.d0 1065 cof(i,j,k,6) = cof(i,j,k,6)+c5 1066 cof(i,j,k,7) = cof(i,j,k,7)+dlz2*c5*alfa 1067 end do 1068 end do 1069 end if 1070 if (nzf.eq.2) then 1071 kbdy = 6 1072 z = zf 1073 k = nz 1074 do j=jst,jfn 1075 y = yc+(j-1)*dly 1076 do i=ist,ifn 1077 x = xa+(i-1)*dlx 1078 call bndyc(kbdy,x,y,alfa,gbdy) 1079 alfmax = max(abs(alfa),alfmax) 1080 c6 = cof(i,j,k,6) 1081 cof(i,j,k,5) = cof(i,j,k,5)+c6 1082 cof(i,j,k,6) = 0.d0 1083 cof(i,j,k,7) = cof(i,j,k,7)-dlz2*c6*alfa 1084 end do 1085 end do 1086 end if 1087c 1088c flag continuous singular elliptic pde if detected 1089c a fatal error for muh3 1090c 1091 if (cemax.eq.0.d0.and.alfmax.eq.0.0) then 1092 if (nxa.eq.0.or.(nxa.eq.2.and.nxb.eq.2)) then 1093 if (nyc.eq.0.or.(nyc.eq.2.and.nyd.eq.2)) then 1094 if (nze.eq.0.or.(nze.eq.2.and.nzf.eq.2)) then 1095 ier =14 1096 end if 1097 end if 1098 end if 1099 end if 1100c 1101c reset cof for specified b.c. 1102c 1103 if (nxa.eq.1) then 1104 i = 1 1105 do j=1,ny 1106 do k=1,nz 1107 do l=1,7 1108 cof(i,j,k,l) = 0.d0 1109 end do 1110 cof(i,j,k,7) = 1.0 1111 end do 1112 end do 1113 end if 1114 if (nxb.eq.1) then 1115 i = nx 1116 do k=1,nz 1117 do j=1,ny 1118 do l=1,7 1119 cof(i,j,k,l) = 0.d0 1120 end do 1121 cof(i,j,k,7) = 1.0 1122 end do 1123 end do 1124 end if 1125 if (nyc.eq.1) then 1126 j = 1 1127 do k=1,nz 1128 do i=1,nx 1129 do l=1,7 1130 cof(i,j,k,l) = 0.d0 1131 end do 1132 cof(i,j,k,7) = 1.0 1133 end do 1134 end do 1135 end if 1136 if (nyd.eq.1) then 1137 j = ny 1138 do i=1,nx 1139 do k=1,nz 1140 do l=1,7 1141 cof(i,j,k,l) = 0.d0 1142 end do 1143 cof(i,j,k,7) = 1.0 1144 end do 1145 end do 1146 end if 1147 if (nze.eq.1) then 1148 k = 1 1149 do j=1,ny 1150 do i=1,nx 1151 do l=1,7 1152 cof(i,j,k,l) = 0.d0 1153 end do 1154 cof(i,j,k,7) = 1.0 1155 end do 1156 end do 1157 end if 1158 if (nzf.eq.1) then 1159 k = nz 1160 do j=1,ny 1161 do i=1,nx 1162 do l=1,7 1163 cof(i,j,k,l) = 0.d0 1164 end do 1165 cof(i,j,k,7) = 1.0 1166 end do 1167 end do 1168 end if 1169c if (klevel.eq.1 .and. method.lt.8) then 1170 if (klevel.eq.1) then 1171 1172c 1173c set coarse grid resolution 1174c 1175 nx = ixp+1 1176 ny = jyq+1 1177 nz = kzr+1 1178 nxny = nx*ny 1179 if (nze .ne. 0) then 1180c 1181c factor non-periodic block matrix 1182c 1183 call lud3(nxny,nx,ny,nz,cof,wk(ibeta),wk(ialfa),iwk,nxa,nyc) 1184 return 1185 else 1186c 1187c factor periodic block matrix 1188c 1189 do k=1,nz-1 1190 call setbet3(nxny,nx,ny,nz,cof,wk(ibeta),k,nxa,nyc) 1191 call setalf3(nxny,nx,ny,nz,cof,wk(ialfa),k) 1192 end do 1193 call lud3p(nxny,nx,ny,nz,cof,wk(ibeta),wk(ialfa),wk(izmat), 1194 + wk(idmat),iwk,nxa,nyc) 1195 return 1196 end if 1197 end if 1198 if (method*(method-8)*(method-9)*(method-10) .eq. 0) return 1199c 1200c set,factor tridiagonal matrices for line relaxation 1201c 1202 if ((method-1)*(method-4)*(method-5)*(method-7).eq.0) then 1203c 1204c line relaxation in x used 1205c 1206 if (nxa.ne.0) then 1207c 1208c set non-periodic tridiagonal matrices in tx and factor 1209c 1210 do i=1,nx 1211 im1 = max0(i-1,1) 1212 do k=1,nz 1213 do j=1,ny 1214 tx(im1,j,k,1) = cof(i,j,k,1) 1215 tx(i,j,k,2) = cof(i,j,k,7) 1216 tx(i,j,k,3) = cof(i,j,k,2) 1217 end do 1218 end do 1219 end do 1220 nynz = ny*nz 1221 call factri(nynz,nx,tx(1,1,1,1),tx(1,1,1,2),tx(1,1,1,3)) 1222 else 1223 if (nx .gt. 3) then 1224c 1225c set "periodic" tridiagonal matrices in tx and factor when nx > 3 1226c 1227 do k=1,nz 1228 do j=1,ny 1229 do i=1,nx-1 1230 tx(i,j,k,1) = cof(i,j,k,1) 1231 tx(i,j,k,2) = cof(i,j,k,7) 1232 tx(i,j,k,3) = cof(i,j,k,2) 1233 end do 1234 end do 1235 end do 1236 nynz = ny*nz 1237 call factrp(nynz,nx,tx,tx(1,1,1,2),tx(1,1,1,3),tx(1,1,1,4), 1238 + tx(1,1,1,5),wk(kps)) 1239 end if 1240 end if 1241 end if 1242 if ((method-2)*(method-4)*(method-6)*(method-7).eq.0) then 1243c 1244c line relaxation in y used 1245c 1246 if (nyc.ne.0) then 1247c 1248c set non-periodic tridiagonal matrices and factor 1249c 1250 do j=1,ny 1251 jm1 = max0(j-1,1) 1252 do k=1,nz 1253 do i=1,nx 1254 ty(jm1,i,k,1) = cof(i,j,k,3) 1255 ty(j,i,k,2) = cof(i,j,k,7) 1256 ty(j,i,k,3) = cof(i,j,k,4) 1257 end do 1258 end do 1259 end do 1260 nxnz = nx*nz 1261 call factri(nxnz,ny,ty(1,1,1,1),ty(1,1,1,2),ty(1,1,1,3)) 1262 else 1263 if (ny .gt. 3) then 1264c 1265c set and factor periodic "tridiagonal" matrices when ny > 3 1266c 1267 do k=1,nz 1268 do i=1,nx 1269 do j=1,ny-1 1270 ty(j,i,k,1) = cof(i,j,k,3) 1271 ty(j,i,k,2) = cof(i,j,k,7) 1272 ty(j,i,k,3) = cof(i,j,k,4) 1273 end do 1274 end do 1275 end do 1276 nxnz = nx*nz 1277 call factrp(nxnz,ny,ty,ty(1,1,1,2),ty(1,1,1,3),ty(1,1,1,4), 1278 + ty(1,1,1,5),wk(kps)) 1279 end if 1280 end if 1281 end if 1282 if ((method-3)*(method-5)*(method-6)*(method-7).eq.0) then 1283c 1284c line relaxation in z used 1285c 1286 if (nze.ne.0) then 1287c 1288c set and factor non-periodic tridiagonal matrices 1289c 1290 do k=1,nz 1291 km1 = max0(k-1,1) 1292 do j=1,ny 1293 do i=1,nx 1294 tz(km1,i,j,1) = cof(i,j,k,5) 1295 tz(k,i,j,2) = cof(i,j,k,7) 1296 tz(k,i,j,3) = cof(i,j,k,6) 1297 end do 1298 end do 1299 end do 1300 nxny = nx*ny 1301 call factri(nxny,nz,tz(1,1,1,1),tz(1,1,1,2),tz(1,1,1,3)) 1302 else 1303 if (nz .gt. 3) then 1304c 1305c set and factor periodic "tridiagonal matrices when nz > 3" 1306c 1307 do j=1,ny 1308 do i=1,nx 1309 do k=1,nz-1 1310 tz(k,i,j,1) = cof(i,j,k,5) 1311 tz(k,i,j,2) = cof(i,j,k,7) 1312 tz(k,i,j,3) = cof(i,j,k,6) 1313 end do 1314 end do 1315 end do 1316 nxny = nx*ny 1317 call factrp(nxny,nz,tz(1,1,1,1),tz(1,1,1,2),tz(1,1,1,3), 1318 + tz(1,1,1,4),tz(1,1,1,5),wk(kps)) 1319 end if 1320 end if 1321 end if 1322 return 1323 end 1324 1325 subroutine adjmh3(nx,ny,nz,phi,cof,bndyc,coef) 1326c 1327c adjust rhs for solution in cof(i,j,k,8) on non-initial calls 1328c (i.e., values in cof have not changed) 1329c 1330 implicit none 1331 integer nx,ny,nz 1332 double precision phi(0:nx+1,0:ny+1,0:nz+1),cof(nx,ny,nz,8) 1333 integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 1334 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 1335 +kcycle,iprer,ipost,intpol 1336 double precision xa,xb,yc,yd,ze,zf,tolmax,relmax 1337 integer kpbgn,kcbgn,ktxbgn,ktybgn,ktzbgn,nxk,nyk,nzk,ngrid, 1338 + klevel,kcur,kps 1339 common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 1340 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 1341 +kcycle,iprer,ipost,intpol 1342 common/fmud3/xa,xb,yc,yd,ze,zf,tolmax,relmax 1343 common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), 1344 +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps 1345 integer ibeta,ialfa,izmat,idmat 1346 common/muh3c/ibeta,ialfa,izmat,idmat 1347 double precision dlx,dly,dlz,dlx2,dly2,dlz2,dlxx,dlyy,dlzz 1348 double precision cxx,cyy,czz,cx,cy,cz,ce,alfa,gbdy,x,y,z,c1,c2,c3, 1349 +c4,c5,c6 1350 integer i,j,k,ist,ifn,jst,jfn,kst,kfn,kbdy 1351 external bndyc,coef 1352c 1353c set current grid increments 1354c 1355 dlx = (xb-xa)/(nx-1) 1356 dlx2 = dlx+dlx 1357 dlxx = dlx*dlx 1358 dly = (yd-yc)/(ny-1) 1359 dly2 = dly+dly 1360 dlyy = dly*dly 1361 dlz = (zf-ze)/(nz-1) 1362 dlz2 = dlz+dlz 1363 dlzz = dlz*dlz 1364c 1365c set x,y,z subscript limits for calls to coef,bndyc 1366c 1367 jst = 1 1368 jfn = ny 1369 ist = 1 1370 ifn = nx 1371 kst = 1 1372 kfn = nz 1373 if (nxa.eq.1) ist = 2 1374 if (nxb.eq.1) ifn = nx-1 1375 if (nyc.eq.1) jst = 2 1376 if (nyd.eq.1) jfn = ny-1 1377 if (nze.eq.1) kst = 2 1378 if (nzf.eq.1) kfn = nz-1 1379c 1380c adjust for derivative b.c. 1381c 1382 if (nxa.eq.2) then 1383 kbdy=1 1384 x=xa 1385 i = 1 1386 do k=kst,kfn 1387 z=ze+(k-1)*dlz 1388 do j=jst,jfn 1389 y=yc+(j-1)*dly 1390 call bndyc(kbdy,y,z,alfa,gbdy) 1391 call coef(x,y,z,cxx,cyy,czz,cx,cy,cz,ce) 1392 cxx = max(cxx,abs(cx)*dlx*0.5d0) 1393 c1 = cxx/dlxx-cx/dlx2 1394 cof(i,j,k,8) = cof(i,j,k,8)+dlx2*c1*gbdy 1395 end do 1396 end do 1397 end if 1398 if (nxb.eq.2) then 1399 kbdy=2 1400 x=xb 1401 i = nx 1402 do k=kst,kfn 1403 z=ze+(k-1)*dlz 1404 do j=jst,jfn 1405 y=yc+(j-1)*dly 1406 call bndyc(kbdy,y,z,alfa,gbdy) 1407 call coef(x,y,z,cxx,cyy,czz,cx,cy,cz,ce) 1408 cxx = max(cxx,abs(cx)*dlx*0.5d0) 1409 c2 = cxx/dlxx+cx/dlx2 1410 cof(i,j,k,8) = cof(i,j,k,8)-dlx2*c2*gbdy 1411 end do 1412 end do 1413 end if 1414 if (nyc.eq.2) then 1415 kbdy = 3 1416 y=yc 1417 j = 1 1418 do k=kst,kfn 1419 z=ze+(k-1)*dlz 1420 do i=ist,ifn 1421 x=xa+(i-1)*dlx 1422 call bndyc(kbdy,x,z,alfa,gbdy) 1423 call coef(x,y,z,cxx,cyy,czz,cx,cy,cz,ce) 1424 cyy = max(cyy,abs(cy)*dly*0.5d0) 1425 c3 = cyy/dlyy-cy/dly2 1426 cof(i,j,k,8) = cof(i,j,k,8)+dly2*c3*gbdy 1427 end do 1428 end do 1429 end if 1430 if (nyd.eq.2) then 1431 kbdy = 4 1432 y=yd 1433 j = ny 1434 do k=kst,kfn 1435 z=ze+(k-1)*dlz 1436 do i=ist,ifn 1437 x=xa+(i-1)*dlx 1438 call bndyc(kbdy,x,z,alfa,gbdy) 1439 call coef(x,y,z,cxx,cyy,czz,cx,cy,cz,ce) 1440 cyy = max(cyy,abs(cy)*dly*0.5d0) 1441 c4 = cyy/dlyy+cy/dly2 1442 cof(i,j,k,8) = cof(i,j,k,8)-dly2*c4*gbdy 1443 end do 1444 end do 1445 end if 1446 if (nze.eq.2) then 1447 kbdy = 5 1448 k = 1 1449 z=ze 1450 do j=jst,jfn 1451 y=yc+(j-1)*dly 1452 do i=ist,ifn 1453 x=xa+(i-1)*dlx 1454 call bndyc(kbdy,x,y,alfa,gbdy) 1455 call coef(x,y,z,cxx,cyy,czz,cx,cy,cz,ce) 1456 czz = max(czz,abs(cz)*dlz*0.5d0) 1457 c5 = czz/dlzz-cz/dlz2 1458 cof(i,j,k,8) = cof(i,j,k,8)+dlz2*c5*gbdy 1459 end do 1460 end do 1461 end if 1462 if (nzf.eq.2) then 1463 kbdy = 6 1464 z=zf 1465 k = nz 1466 do j=jst,jfn 1467 y=yc+(j-1)*dly 1468 do i=ist,ifn 1469 x=xa+(i-1)*dlx 1470 call bndyc(kbdy,x,y,alfa,gbdy) 1471 call coef(x,y,z,cxx,cyy,czz,cx,cy,cz,ce) 1472 czz = max(czz,abs(cz)*dlz*0.5d0) 1473 c6 = czz/dlzz+cz/dlz2 1474 cof(i,j,k,8) = cof(i,j,k,8)-dlz2*c6*gbdy 1475 end do 1476 end do 1477 end if 1478c 1479c set specified b.c. 1480c 1481 if (nxa.eq.1) then 1482 i = 1 1483 do j=1,ny 1484 do k=1,nz 1485 cof(i,j,k,8) = phi(i,j,k) 1486 end do 1487 end do 1488 end if 1489 if (nxb.eq.1) then 1490 i = nx 1491 do j=1,ny 1492 do k=1,nz 1493 cof(i,j,k,8) = phi(i,j,k) 1494 end do 1495 end do 1496 end if 1497 if (nyc.eq.1) then 1498 j = 1 1499 do k=1,nz 1500 do i=1,nx 1501 cof(i,j,k,8) = phi(i,j,k) 1502 end do 1503 end do 1504 end if 1505 if (nyd.eq.1) then 1506 j = ny 1507 do k=1,nz 1508 do i=1,nx 1509 cof(i,j,k,8) = phi(i,j,k) 1510 end do 1511 end do 1512 end if 1513 if (nze.eq.1) then 1514 k = 1 1515 do j=1,ny 1516 do i=1,nx 1517 cof(i,j,k,8) = phi(i,j,k) 1518 end do 1519 end do 1520 end if 1521 if (nzf.eq.1) then 1522 k = nz 1523 do j=1,ny 1524 do i=1,nx 1525 cof(i,j,k,8) = phi(i,j,k) 1526 end do 1527 end do 1528 end if 1529 return 1530 end 1531 1532 subroutine relmh3(wk) 1533c 1534c use point or line relaxation in the x and/or y and/or z 1535c or planar relaxation in the x,y or x,z or y,z planes 1536c 1537 implicit none 1538 double precision wk(*) 1539 integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 1540 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 1541 +kcycle,iprer,ipost,intpol 1542 double precision xa,xb,yc,yd,ze,zf,tolmax,relmax 1543 integer kpbgn,kcbgn,ktxbgn,ktybgn,ktzbgn,nxk,nyk,nzk,ngrid, 1544 + klevel,kcur,kps 1545 common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 1546 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 1547 +kcycle,iprer,ipost,intpol 1548 common/fmud3/xa,xb,yc,yd,ze,zf,tolmax,relmax 1549 common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), 1550 +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps 1551 integer ibeta,ialfa,izmat,idmat 1552 common/muh3c/ibeta,ialfa,izmat,idmat 1553 integer nx,ny,nz,ip,ic,m,itx,ity,itz 1554 nx = nxk(klevel) 1555 ny = nyk(klevel) 1556 nz = nzk(klevel) 1557 ip = kpbgn(klevel) 1558 ic = kcbgn(klevel) 1559 if (method.eq.0) then 1560c 1561c gauss-seidel pointwise red/black relaxation 1562c 1563 call relmh3p(nx,ny,nz,wk(ip),wk(ic)) 1564 return 1565 end if 1566 itx = ktxbgn(klevel) 1567 ity = ktybgn(klevel) 1568 itz = ktzbgn(klevel) 1569 m = method 1570c 1571c check for line relaxation(s) (in combinations) 1572c 1573 if ((m-1)*(m-4)*(m-5)*(m-7) .eq. 0 ) then 1574c 1575c line - x relaxation 1576c 1577 if (nxa .ne. 0 .or. nx .gt. 3) then 1578 itx = ktxbgn(klevel) 1579 call slxmd3(nx,ny,nz,wk(ip),wk(ic),wk(itx),wk(kps),nxa,nyc,nze) 1580 else 1581c 1582c replace by point if x-periodic and nx=3 1583c 1584 call relmh3p(nx,ny,nz,wk(ip),wk(ic)) 1585 end if 1586 if (method .eq. 1) return 1587 end if 1588 1589 if ((m-2)*(m-4)*(m-6)*(m-7) .eq. 0 ) then 1590c 1591c line - y relaxation 1592c 1593 if (nyc .ne. 0 .or. ny .gt. 3) then 1594 ity = ktybgn(klevel) 1595 call slymd3(nx,ny,nz,wk(ip),wk(ic),wk(ity),wk(kps),nxa,nyc,nze) 1596 else 1597c 1598c replace by point if y-periodic and ny=3 1599c 1600 call relmh3p(nx,ny,nz,wk(ip),wk(ic)) 1601 end if 1602 if ((m-2)*(m-4) .eq. 0) return 1603 end if 1604 1605 if ((m-3)*(m-5)*(m-6)*(m-7) .eq. 0 ) then 1606c 1607c line - z relaxation 1608c 1609 if (nze .ne. 0 .or. nz .gt. 3) then 1610 itz = ktzbgn(klevel) 1611 call slzmd3(nx,ny,nz,wk(ip),wk(ic),wk(itz),wk(kps),nxa,nyc,nze) 1612 else 1613c 1614c replace by point if z-periodic and nz=3 1615c 1616 call relmh3p(nx,ny,nz,wk(ip),wk(ic)) 1617 end if 1618 return 1619 end if 1620c 1621c planar relaxation 1622c 1623 if (method.eq.8) then 1624c 1625c planar relaxation in the x,y plane 1626c 1627 call planxy(wk) 1628 else if (method.eq.9) then 1629c 1630c planar relaxation in the x,z plane 1631c 1632 call planxz(wk) 1633 else if (method.eq.10) then 1634c 1635c planar relaxation in the y,z plane 1636c 1637 call planyz(wk) 1638 end if 1639 return 1640 end 1641 1642 subroutine relmh3p(nx,ny,nz,phi,cof) 1643c 1644c gauss-seidel point relaxation with red/black ordering 1645c in three dimensions for nonseparable pde 1646c 1647 implicit none 1648 integer nx,ny,nz,i,j,k,nper 1649 integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 1650 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 1651 +kcycle,iprer,ipost,intpol 1652 common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 1653 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 1654 +kcycle,iprer,ipost,intpol 1655 double precision phi(0:nx+1,0:ny+1,0:nz+1),cof(nx,ny,nz,8) 1656c 1657c set periodic b.c. indicator 1658c 1659 nper = nxa*nyc*nze 1660c 1661c set periodic boundaries as necessary 1662c 1663 if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze) 1664c 1665c relax in order: 1666c (1) red (x,y) on odd z planes 1667c (2) black (x,y) on even z planes 1668c (3) black (x,y) on odd z planes 1669c (4) red (x,y) on even z planes 1670c 1671!$OMP PARALLEL DO PRIVATE(i,j,k), SHARED(phi,cof,nx,ny,nz) 1672 do k=1,nz,2 1673c 1674c red (x,y) points on odd z planes 1675c 1676 do i=1,nx,2 1677 do j=1,ny,2 1678 phi(i,j,k) = (cof(i,j,k,8) - ( 1679 + cof(i,j,k,1)*phi(i-1,j,k)+ 1680 + cof(i,j,k,2)*phi(i+1,j,k)+ 1681 + cof(i,j,k,3)*phi(i,j-1,k)+ 1682 + cof(i,j,k,4)*phi(i,j+1,k)+ 1683 + cof(i,j,k,5)*phi(i,j,k-1)+ 1684 + cof(i,j,k,6)*phi(i,j,k+1))) 1685 + /cof(i,j,k,7) 1686 end do 1687 end do 1688 do i=2,nx,2 1689 do j=2,ny,2 1690 phi(i,j,k) = (cof(i,j,k,8) - ( 1691 + cof(i,j,k,1)*phi(i-1,j,k)+ 1692 + cof(i,j,k,2)*phi(i+1,j,k)+ 1693 + cof(i,j,k,3)*phi(i,j-1,k)+ 1694 + cof(i,j,k,4)*phi(i,j+1,k)+ 1695 + cof(i,j,k,5)*phi(i,j,k-1)+ 1696 + cof(i,j,k,6)*phi(i,j,k+1))) 1697 + /cof(i,j,k,7) 1698 end do 1699 end do 1700 end do 1701!$OMP END PARALLEL DO 1702 if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze) 1703c 1704c black (x,y) points on even z planes 1705c 1706!$OMP PARALLEL DO PRIVATE(i,j,k), SHARED(phi,cof,nx,ny,nz) 1707 do k=2,nz,2 1708 do i=1,nx,2 1709 do j=2,ny,2 1710 phi(i,j,k) = (cof(i,j,k,8) - ( 1711 + cof(i,j,k,1)*phi(i-1,j,k)+ 1712 + cof(i,j,k,2)*phi(i+1,j,k)+ 1713 + cof(i,j,k,3)*phi(i,j-1,k)+ 1714 + cof(i,j,k,4)*phi(i,j+1,k)+ 1715 + cof(i,j,k,5)*phi(i,j,k-1)+ 1716 + cof(i,j,k,6)*phi(i,j,k+1))) 1717 + /cof(i,j,k,7) 1718 end do 1719 end do 1720 do i=2,nx,2 1721 do j=1,ny,2 1722 phi(i,j,k) = (cof(i,j,k,8) - ( 1723 + cof(i,j,k,1)*phi(i-1,j,k)+ 1724 + cof(i,j,k,2)*phi(i+1,j,k)+ 1725 + cof(i,j,k,3)*phi(i,j-1,k)+ 1726 + cof(i,j,k,4)*phi(i,j+1,k)+ 1727 + cof(i,j,k,5)*phi(i,j,k-1)+ 1728 + cof(i,j,k,6)*phi(i,j,k+1))) 1729 + /cof(i,j,k,7) 1730 end do 1731 end do 1732 end do 1733!$OMP END PARALLEL DO 1734 if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze) 1735c 1736c black (x,y) points on odd z planes 1737c 1738!$OMP PARALLEL DO PRIVATE(i,j,k), SHARED(phi,cof,nx,ny,nz) 1739 do k=1,nz,2 1740 do i=1,nx,2 1741 do j=2,ny,2 1742 phi(i,j,k) = (cof(i,j,k,8) - ( 1743 + cof(i,j,k,1)*phi(i-1,j,k)+ 1744 + cof(i,j,k,2)*phi(i+1,j,k)+ 1745 + cof(i,j,k,3)*phi(i,j-1,k)+ 1746 + cof(i,j,k,4)*phi(i,j+1,k)+ 1747 + cof(i,j,k,5)*phi(i,j,k-1)+ 1748 + cof(i,j,k,6)*phi(i,j,k+1))) 1749 + /cof(i,j,k,7) 1750 end do 1751 end do 1752 do i=2,nx,2 1753 do j=1,ny,2 1754 phi(i,j,k) = (cof(i,j,k,8) - ( 1755 + cof(i,j,k,1)*phi(i-1,j,k)+ 1756 + cof(i,j,k,2)*phi(i+1,j,k)+ 1757 + cof(i,j,k,3)*phi(i,j-1,k)+ 1758 + cof(i,j,k,4)*phi(i,j+1,k)+ 1759 + cof(i,j,k,5)*phi(i,j,k-1)+ 1760 + cof(i,j,k,6)*phi(i,j,k+1))) 1761 + /cof(i,j,k,7) 1762 end do 1763 end do 1764 end do 1765!$OMP END PARALLEL DO 1766 if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze) 1767c 1768c red (x,y) points on even z planes 1769c 1770!$OMP PARALLEL DO PRIVATE(i,j,k), SHARED(phi,cof,nx,ny,nz) 1771 do k=2,nz,2 1772 do i=1,nx,2 1773 do j=1,ny,2 1774 phi(i,j,k) = (cof(i,j,k,8) - ( 1775 + cof(i,j,k,1)*phi(i-1,j,k)+ 1776 + cof(i,j,k,2)*phi(i+1,j,k)+ 1777 + cof(i,j,k,3)*phi(i,j-1,k)+ 1778 + cof(i,j,k,4)*phi(i,j+1,k)+ 1779 + cof(i,j,k,5)*phi(i,j,k-1)+ 1780 + cof(i,j,k,6)*phi(i,j,k+1))) 1781 + /cof(i,j,k,7) 1782 end do 1783 end do 1784 do i=2,nx,2 1785 do j=2,ny,2 1786 phi(i,j,k) = (cof(i,j,k,8) - ( 1787 + cof(i,j,k,1)*phi(i-1,j,k)+ 1788 + cof(i,j,k,2)*phi(i+1,j,k)+ 1789 + cof(i,j,k,3)*phi(i,j-1,k)+ 1790 + cof(i,j,k,4)*phi(i,j+1,k)+ 1791 + cof(i,j,k,5)*phi(i,j,k-1)+ 1792 + cof(i,j,k,6)*phi(i,j,k+1))) 1793 + /cof(i,j,k,7) 1794 end do 1795 end do 1796 end do 1797!$OMP END PARALLEL DO 1798c 1799c final set of periodic virtual boundaries if necessary 1800c 1801 if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze) 1802 return 1803 end 1804 1805 subroutine lud3(nxny,nx,ny,nz,cof,beta,alfa,index,nxa,nyc) 1806c 1807c decompose nonperiodic block coefficient matrix 1808c coming from 3-d PDE discretization 1809c 1810 implicit none 1811 integer nxny,nx,ny,nz,nxa,nyc 1812 double precision cof(nx,ny,nz,8),beta(nxny,nxny,*) 1813 double precision alfa(nxny,nxny,*) 1814 integer index(nxny,nz) 1815 integer iz,i1,kcur,ii,ll,il,jl,km1 1816 iz = 0 1817 i1 = 1 1818c 1819c set and factor umat(1) in beta(1) 1820c 1821 kcur = 1 1822 call setbet3(nxny,nx,ny,nz,cof,beta,kcur,nxa,nyc) 1823 call sgfa(beta,nxny,nxny,index,iz) 1824 do kcur=2,nz 1825 call setalf3(nxny,nx,ny,nz,cof,alfa,kcur) 1826 call transp(nxny,alfa(1,1,kcur)) 1827 km1 = kcur-1 1828 do ll=1,nxny 1829 call sgsl(beta(1,1,km1),nxny,nxny,index(1,km1), 1830 + alfa(1,ll,kcur),i1) 1831 end do 1832 call transp(nxny,alfa(1,1,kcur)) 1833c 1834C set beta(kcur)=beta(kcur)-alfa(kcur)*gama(kcur-1) 1835c 1836 call setbet3(nxny,nx,ny,nz,cof,beta,kcur,nxa,nyc) 1837 do ii=1,nxny 1838 do ll=1,nxny 1839 jl = ((ll-1)/nx) + 1 1840 il = ll - (jl-1)*nx 1841 beta(ii,ll,kcur)=beta(ii,ll,kcur)-alfa(ii,ll,kcur)* 1842 + cof(il,jl,kcur-1,6) 1843 end do 1844 end do 1845c 1846c factor current beta for next pass 1847c 1848 iz = 0 1849 call sgfa(beta(1,1,kcur),nxny,nxny,index(1,kcur),iz) 1850 end do 1851c 1852C matrix factorization now complete 1853c 1854 return 1855 end 1856 1857 subroutine dir3(nxny,nx,ny,nz,phi,cof,beta,alfa,phxy,index) 1858c 1859c direct solve at coarsest grid 1860c 1861 implicit none 1862 integer nxny,nx,ny,nz,index(nxny,ny) 1863 double precision phi(0:nx+1,0:ny+1,0:nz+1),cof(nx,ny,nz,8) 1864 double precision phxy(nxny) 1865 double precision beta(nxny,nxny,*),alfa(nxny,nxny,*) 1866c forward sweep 1867 call for3(nxny,nx,ny,nz,phi,cof(1,1,1,8),alfa) 1868c backward sweep 1869 call bkw3(nxny,nx,ny,nz,phi,cof,beta,phxy,index) 1870 return 1871 end 1872 1873 subroutine for3(nxny,nx,ny,nz,phi,frhs,alfa) 1874c 1875c forward sweep 1876c 1877 implicit none 1878 integer nxny,nx,ny,nz,i,j,k,ii,ll,il,jl 1879 double precision phi(0:nx+1,0:ny+1,0:nz+1) 1880 double precision frhs(nx,ny,nz),alfa(nxny,nxny,*) 1881 double precision sum 1882 do k=1,nz 1883 do j=1,ny 1884 do i=1,nx 1885 phi(i,j,k) = frhs(i,j,k) 1886 end do 1887 end do 1888 end do 1889 do k=2,nz 1890c 1891C SET PHI(k)=PHI(k)-alfa(k)*PHI(k-1) 1892c 1893 do ii=1,nxny 1894 j = ((ii-1)/nx)+1 1895 i = ii - (j-1)*nx 1896 sum=0.d0 1897 do ll=1,nxny 1898 jl = ((ll-1)/nx)+1 1899 il = ll - (jl-1)*nx 1900 sum = sum+alfa(ii,ll,k)*phi(il,jl,k-1) 1901 end do 1902 phi(i,j,k) = phi(i,j,k)-sum 1903 end do 1904 end do 1905 return 1906 end 1907 1908 subroutine bkw3(nxny,nx,ny,nz,phi,cof,beta,phxy,index) 1909C 1910C backward sweep 1911c 1912 implicit none 1913 integer nxny,nx,ny,nz,index(nxny,nz) 1914 double precision phi(0:nx+1,0:ny+1,0:nz+1),cof(nx,ny,nz,8) 1915 double precision beta(nxny,nxny,*),phxy(nxny) 1916 integer iz,i,j,k,kb,ij 1917C 1918C solve beta*phi(nz)=phi(nz) 1919C 1920 iz = 0 1921 do j=1,ny 1922 do i=1,nx 1923 ij = (j-1)*nx+i 1924 phxy(ij) = phi(i,j,nz) 1925 end do 1926 end do 1927 call sgsl(beta(1,1,nz),nxny,nxny,index(1,nz),phxy,iz) 1928 do j=1,ny 1929 do i=1,nx 1930 ij = (j-1)*nx+i 1931 phi(i,j,nz) = phxy(ij) 1932 end do 1933 end do 1934 do kb=2,nz 1935 k = nz-kb+1 1936C 1937C set current rhs 1938C 1939 do j=1,ny 1940 do i=1,nx 1941 phi(i,j,k) = phi(i,j,k) - cof(i,j,k,6)*phi(i,j,k+1) 1942 end do 1943 end do 1944C 1945C solve beta(k)*phi(k)=phi(k) 1946C 1947 do j=1,ny 1948 do i=1,nx 1949 ij = (j-1)*nx+i 1950 phxy(ij) = phi(i,j,k) 1951 end do 1952 end do 1953 call sgsl(beta(1,1,k),nxny,nxny,index(1,k),phxy,iz) 1954 do j=1,ny 1955 do i=1,nx 1956 ij = (j-1)*nx+i 1957 phi(i,j,k) = phxy(ij) 1958 end do 1959 end do 1960 end do 1961 return 1962 end 1963 1964 subroutine lud3p(nxny,nx,ny,nz,cof,beta,alfa,zmat,dmat,index, 1965 + nxa,nyc) 1966C 1967C DO LUD DECOMPOSITION OF BLOCK periodic TRIDIAGONAL MATRIX EQUATION 1968c 1969 implicit none 1970 integer nxny,nx,ny,nz,index(nxny,nz),nxa,nyc 1971 double precision cof(nx,ny,nz,8),alfa(nxny,nxny,*), 1972 +beta(nxny,nxny,*) 1973 double precision dmat(nxny,nxny,*),zmat(nxny,nxny,*),sum 1974 integer iz,kcur,ii,jj,ll,i1,km1,il,jl,i,j,ij 1975 kcur = 1 1976c 1977c set dmat(1)=alfa(1) 1978c 1979 call setalf3(nxny,nx,ny,nz,cof,alfa,kcur) 1980 do ii=1,nxny 1981 do ll=1,nxny 1982 dmat(ii,ll,1) = alfa(ii,ll,kcur) 1983 end do 1984 end do 1985c 1986c factor umat(1) in beta(1) 1987c 1988 call setbet3(nxny,nx,ny,nz,cof,beta,kcur,nxa,nyc) 1989 iz = 0 1990 call sgfa(beta(1,1,1),nxny,nxny,index(1,1),iz) 1991 do kcur=2,nz-2 1992c 1993c solve transpose of lmat(kcur)umat(kcur-1)=alfa(kcur) in alfa(kcur) 1994c 1995 call setalf3(nxny,nx,ny,nz,cof,alfa,kcur) 1996c 1997c transpose alfa 1998c 1999 call transp(nxny,alfa(1,1,kcur)) 2000c 2001c solve transpose equation 2002c 2003 km1 = kcur-1 2004 i1 = 1 2005 do ll=1,nxny 2006 call sgsl(beta(1,1,km1),nxny,nxny,index(1,km1),alfa(1,ll,kcur) 2007 + ,i1) 2008 end do 2009c 2010c re-transpose solution in alfa 2011c 2012 call transp(nxny,alfa(1,1,kcur)) 2013c 2014C set umat(kcur) = beta(kcur) - alfa(kcur)*gama(kcur-1) in beta(kcur) 2015c 2016 call setbet3(nxny,nx,ny,nz,cof,beta,kcur,nxa,nyc) 2017 do ii=1,nxny 2018 do ll=1,nxny 2019 jl = ((ll-1)/nx)+1 2020 il = ll-(jl-1)*nx 2021 beta(ii,ll,kcur)=beta(ii,ll,kcur)-alfa(ii,ll,kcur)* 2022 + cof(il,jl,kcur-1,6) 2023 end do 2024 end do 2025c 2026c factor current beta(1,1,kcur) for next pass 2027c 2028 call sgfa(beta(1,1,kcur),nxny,nxny,index(1,kcur),iz) 2029c 2030c set dmat(kcur) = -alfa(kcur)*dmat(kcur-1) 2031c 2032 do ii=1,nxny 2033 do jj=1,nxny 2034 dmat(ii,jj,kcur) = 0.d0 2035 do ll=1,nxny 2036 dmat(ii,jj,kcur)=dmat(ii,jj,kcur)-alfa(ii,ll,kcur)* 2037 + dmat(ll,jj,kcur-1) 2038 end do 2039 end do 2040 end do 2041 if (kcur .eq. nz-2) then 2042c 2043c adjust dmat(nz-2) = gama(nz-2)-alfa(nz-2)*dmat(nz-2) 2044c 2045 do j=1,ny 2046 do i=1,nx 2047 ij = (j-1)*nx+i 2048 dmat(ij,ij,kcur) = cof(i,j,kcur,6)+dmat(ij,ij,kcur) 2049 end do 2050 end do 2051 end if 2052 end do 2053c 2054c final phase with periodic factorization 2055c 2056c solve transpose of zmat(1) beta(1) = gama(nz-1) 2057c 2058 do ii=1,nxny 2059 j = ((ii-1)/nx)+1 2060 i = ii-(j-1)*nx 2061 do jj=1,nxny 2062 zmat(jj,ii,1) = 0.d0 2063 end do 2064 zmat(ii,ii,1) = cof(i,j,nz-1,6) 2065 end do 2066 do ll=1,nxny 2067 call sgsl(beta(1,1,1),nxny,nxny,index(1,1),zmat(1,ll,1),i1) 2068 end do 2069 call transp(nxny,zmat(1,1,1)) 2070c 2071c solve transpose of zmat(k) umat(k) = -zmat(k-1) gama(k-1) 2072c 2073 do kcur = 2,nz-3 2074 do ii=1,nxny 2075 do jj=1,nxny 2076 j = ((jj-1)/nx)+1 2077 i = jj-(j-1)*nx 2078 zmat(jj,ii,kcur) = -zmat(ii,jj,kcur-1)*cof(i,j,kcur-1,6) 2079 end do 2080 end do 2081 do ll=1,nxny 2082 call sgsl(beta(1,1,kcur),nxny,nxny,index(1,kcur), 2083 + zmat(1,ll,kcur),i1) 2084 end do 2085 call transp(nxny,zmat(1,1,kcur)) 2086 end do 2087c 2088c solve transpose of zmat(nz-2)umat(nz-2)=alfa(nz-1)-zmat(nz-3)gama(nz-3) 2089c 2090 call setalf3(nxny,nx,ny,nz,cof,alfa,nz-1) 2091 do ii=1,nxny 2092 do jj=1,nxny 2093 j = ((jj-1)/nx)+1 2094 i = jj-(j-1)*nx 2095 zmat(jj,ii,nz-2) = alfa(ii,jj,nz-1)-zmat(ii,jj,nz-3)* 2096 + cof(i,j,nz-3,6) 2097 end do 2098 end do 2099 do ll=1,nxny 2100 call sgsl(beta(1,1,nz-2),nxny,nxny,index(1,nz-2), 2101 + zmat(1,ll,nz-2),i1) 2102 end do 2103 call transp(nxny,zmat(1,1,nz-2)) 2104c 2105c set umat(nz-1) = beta(nz-1)-(zmat(1)*dmat(1)+...+zmat(nz-2)*dmat(nz-2)) 2106c in beta(nz-1) 2107c 2108 call setbet3(nxny,nx,ny,nz,cof,beta,nz-1,nxa,nyc) 2109 do ii=1,nxny 2110 do jj=1,nxny 2111 sum = 0.d0 2112 do kcur=1,nz-2 2113 do ll=1,nxny 2114 sum = sum + zmat(ii,ll,kcur)*dmat(ll,jj,kcur) 2115 end do 2116 end do 2117 beta(ii,jj,nz-1) = beta(ii,jj,nz-1) - sum 2118 end do 2119 end do 2120c 2121c factor bmat(nz-1) for backward sweep 2122c 2123 call sgfa(beta(1,1,nz-1),nxny,nxny,index(1,nz-1),iz) 2124c 2125c lud is now complete 2126c 2127 return 2128 end 2129 2130 subroutine dir3p(nxny,nx,ny,nz,phi,cof,beta,alfa,zmat,dmat, 2131 + phxy,index) 2132c 2133C direct solution of periodic block matrix at coarsest level 2134c 2135 implicit none 2136 integer nxny,nx,ny,nz,index(nxny,nz) 2137 double precision phi(0:nx+1,0:ny+1,0:nz+1),cof(nx,ny,nz,8) 2138 double precision beta(nxny,nxny,*),alfa(nxny,nxny,*) 2139 double precision zmat(nxny,nxny,*), dmat(nxny,nxny,*) 2140 double precision phxy(nxny) 2141c forward sweep for periodic 2142 call for3p(nxny,nx,ny,nz,phi,cof(1,1,1,8),alfa,zmat) 2143c backward sweep for periodic 2144 call bkw3p(nxny,nx,ny,nz,phi,cof,beta,dmat,phxy,index) 2145 return 2146 end 2147 2148 subroutine for3p(nxny,nx,ny,nz,phi,frhs,alfa,zmat) 2149C 2150C DO FORWARD SWEEP PHASE OF DIRECT METHOD IN SOLUTION OF 2151C periodic FACTORED BLOCK TRIDIAGONAL SYSTEM 2152C 2153 implicit none 2154 integer nxny,nx,ny,nz 2155 double precision phi(0:nx+1,0:ny+1,0:nz+1) 2156 double precision frhs(nx,ny,nz) 2157 double precision alfa(nxny,nxny,*),zmat(nxny,nxny,*),sum 2158 integer i,j,k,kcur,ii,ll,il,jl 2159C SET RHS IN PHI ADJUSTED AT DIRICHLET BOUNDARIES 2160 do k=1,nz-1 2161 do j=1,ny 2162 do i=1,nx 2163 phi(i,j,k) = frhs(i,j,k) 2164 end do 2165 end do 2166 end do 2167 do kcur=2,nz-2 2168 do ii=1,nxny 2169 j = ((ii-1)/nx)+1 2170 i = ii-(j-1)*nx 2171 sum = 0.d0 2172 do ll=1,nxny 2173 jl = ((ll-1)/nx)+1 2174 il = ll - (jl-1)*nx 2175 sum=sum+alfa(ii,ll,kcur)*phi(il,jl,kcur-1) 2176 end do 2177 phi(i,j,kcur) = phi(i,j,kcur)-sum 2178 end do 2179 end do 2180c 2181c solve: 2182c zmat(1)*phi(1)+...+zmat(nz-2)*phi(nz-2) + phi(nz-1) = f(nz-1) 2183c 2184 do ii=1,nxny 2185 j = ((ii-1)/nx)+1 2186 i = ii-(j-1)*nx 2187 sum = 0.d0 2188 do k=1,nz-2 2189 do ll=1,nxny 2190 jl = ((ll-1)/nx)+1 2191 il = ll - (jl-1)*nx 2192 sum = sum + zmat(ii,ll,k)*phi(il,jl,k) 2193 end do 2194 end do 2195 phi(i,j,nz-1) = phi(i,j,nz-1) - sum 2196 end do 2197 return 2198 end 2199 2200 subroutine bkw3p(nxny,nx,ny,nz,phi,cof,beta,dmat,phxy,index) 2201C 2202C DO BACKWARD SWEEP 2203c 2204 implicit none 2205 integer nxny,nx,ny,nz,index(nxny,nz) 2206 double precision phi(0:nx+1,0:ny+1,0:nz+1),cof(nx,ny,nz,8) 2207 double precision beta(nxny,nxny,nz),dmat(nxny,nxny,nz) 2208 double precision phxy(nxny),sum 2209 integer iz,ii,ll,il,jl,k,kb,i,j,ij 2210C 2211C solve beta(nz-1)*phi(nz-1) = phi(nz-1) 2212C 2213 iz = 0 2214 do j=1,ny 2215 do i=1,nx 2216 ij = (j-1)*nx+i 2217 phxy(ij) = phi(i,j,nz-1) 2218 end do 2219 end do 2220 2221 call sgsl(beta(1,1,nz-1),nxny,nxny,index(1,nz-1),phxy,iz) 2222 do j=1,ny 2223 do i=1,nx 2224 ij = (j-1)*nx+i 2225 phi(i,j,nz-1) = phxy(ij) 2226 end do 2227 end do 2228c 2229c solve beta(nz-2)*phi(nz-2) = phi(nz-2)-dmat(nz-2)*phi(nz-1) 2230c 2231 do ii=1,nxny 2232 j = ((ii-1)/nx)+1 2233 i = ii-(j-1)*nx 2234 sum = 0.d0 2235 do ll=1,nxny 2236 jl = ((ll-1)/nx)+1 2237 il = ll - (jl-1)*nx 2238 sum = sum + dmat(ii,ll,nz-2)*phi(il,jl,nz-1) 2239 end do 2240 phi(i,j,nz-2) = phi(i,j,nz-2) - sum 2241 end do 2242 do j=1,ny 2243 do i=1,nx 2244 ij = (j-1)*nx+i 2245 phxy(ij) = phi(i,j,nz-2) 2246 end do 2247 end do 2248 call sgsl(beta(1,1,nz-2),nxny,nxny,index(1,nz-2),phxy,iz) 2249 do j=1,ny 2250 do i=1,nx 2251 ij = (j-1)*nx+i 2252 phi(i,j,nz-2) = phxy(ij) 2253 end do 2254 end do 2255c 2256c solve beta(k)*phi(k) = phi(k) - gama(k)*phi(k+1)-dmat(k)*phi(nz-1) 2257c k=nz-3,...,1 2258c 2259 do kb=4,nz 2260 k = nz-kb+1 2261 do ii=1,nxny 2262 j = ((ii-1)/nx)+1 2263 i = ii-(j-1)*nx 2264 sum = 0.d0 2265 do ll=1,nxny 2266 jl = ((ll-1)/nx)+1 2267 il = ll - (jl-1)*nx 2268 sum = sum+dmat(ii,ll,k)*phi(il,jl,nz-1) 2269 end do 2270 phi(i,j,k) = phi(i,j,k) - sum - cof(i,j,k,6)*phi(i,j,k+1) 2271 end do 2272 do j=1,ny 2273 do i=1,nx 2274 ij = (j-1)*nx+i 2275 phxy(ij) = phi(i,j,k) 2276 end do 2277 end do 2278 call sgsl(beta(1,1,k),nxny,nxny,index(1,k),phxy,iz) 2279 do j=1,ny 2280 do i=1,nx 2281 ij = (j-1)*nx+i 2282 phi(i,j,k) = phxy(ij) 2283 end do 2284 end do 2285 end do 2286c 2287c set k=nz by periodicity 2288c 2289 do j=1,ny 2290 do i=1,nx 2291 phi(i,j,nz) = phi(i,j,1) 2292 end do 2293 end do 2294 return 2295 end 2296 2297 subroutine setbet3(nxny,nx,ny,nz,cof,beta,kcur,nxa,nyc) 2298c 2299c set diagonal matrix on block for k = kcur 2300c 2301 implicit none 2302 integer nxny,nx,ny,nz,nxa,nyc 2303 double precision cof(nx,ny,nz,8),beta(nxny,nxny,*) 2304 integer k,kcur,ii,ll,ij,jj,i,j 2305 k = kcur 2306 do ii = 1,nxny 2307 do ll=1,nxny 2308 beta(ii,ll,k)=0.d0 2309 end do 2310 end do 2311 do j=1,ny 2312 do i=1,nx 2313 ij = (j-1)*nx+i 2314 beta(ij,ij,k) = cof(i,j,k,7) 2315 end do 2316 end do 2317 do i=2,nx 2318 do j=1,ny 2319 ij = (j-1)*nx+i 2320 beta(ij,ij-1,k) = cof(i,j,k,1) 2321 end do 2322 end do 2323 do i=1,nx-1 2324 do j=1,ny 2325 ij = (j-1)*nx+i 2326 beta(ij,ij+1,k) = cof(i,j,k,2) 2327 end do 2328 end do 2329 do j=2,ny 2330 do i=1,nx 2331 ij = (j-1)*nx+i 2332 beta(ij,ij-nx,k) = cof(i,j,k,3) 2333 end do 2334 end do 2335 do j=1,ny-1 2336 do i=1,nx 2337 ij = (j-1)*nx+i 2338 beta(ij,ij+nx,k) = cof(i,j,k,4) 2339 end do 2340 end do 2341c 2342c adjust for periodic b.c in x and/or y as necessary 2343c 2344 if (nxa .eq. 0) then 2345 do j=1,ny 2346 ii = (j-1)*nx+1 2347 jj = (j-1)*nx+nx-1 2348 beta(ii,jj,k) = cof(1,j,k,1) 2349 ii = (j-1)*nx+nx 2350 jj = (j-1)*nx+2 2351 beta(ii,jj,k) = cof(nx,j,k,2) 2352 end do 2353 end if 2354 2355 if (nyc.eq.0) then 2356 do i=1,nx 2357 ii = i 2358 jj = (ny-2)*nx+i 2359 beta(ii,jj,k) = cof(i,1,k,3) 2360 ii = (ny-1)*nx+i 2361 jj = nx + i 2362 beta(ii,jj,k) = cof(i,ny,k,4) 2363 end do 2364 end if 2365 return 2366 end 2367 2368 subroutine setalf3(nxny,nx,ny,nz,cof,alfa,kcur) 2369C 2370C set block subdiagonal matrix for k = kcur 2371C 2372 implicit none 2373 integer nxny,nx,ny,nz,kcur 2374 double precision cof(nx,ny,nz,8),alfa(nxny,nxny,*) 2375 integer k,j,i,ij,ll 2376 k = kcur 2377 do j=1,ny 2378 do i=1,nx 2379 ij = (j-1)*nx+i 2380 do ll=1,nxny 2381 alfa(ij,ll,k)=0.d0 2382 end do 2383 alfa(ij,ij,k)=cof(i,j,k,5) 2384 end do 2385 end do 2386 return 2387 end 2388