1 subroutine Bflux ( y, ac, u, x, 2 & shp, shgl, shpb, 3 & shglb, ilwork, iBC, 4 & BC, iper , wallssVec) 5c 6c---------------------------------------------------------------------- 7c 8c This routine : 9c 1. computes the boundary fluxes 10c 2. prints the results in the file [FLUX.lstep] 11c 12c output: 13c in file flux.<lstep>.n (similar to restart file): 14c machin nshg lstep 15c normal_1 ... normal_nsd ! outward normal direction 16c tau_1n ... tau_nsd n ! boundary viscous flux 17c 18c Zdenek Johan, Summer 1991. 19c---------------------------------------------------------------------- 20c 21 22 use pointer_data 23 24 include "common.h" 25 include "mpif.h" 26 27 character*10 cname2 28 29 real*8 y(nshg,ndof), ac(nshg,ndof), 30 & u(nshg,nsd), x(numnp,nsd) 31 dimension iBC(nshg), 32 & BC(nshg,ndofBC), 33 & iper(nshg) 34 35 real*8 shp(MAXTOP,maxsh,MAXQPT), 36 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 37 & shpb(MAXTOP,maxsh,MAXQPT), 38 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 39c 40c 41 real*8 flxres(nshg,nflow), 42 & flxLHS(nshg,1), flxnrm(nshg,nsd), 43 & temp(nshg), rtmp(nshg,ndof), 44 & flxTot(nflow), wallssVec(nshg,3) 45 46 real*8 qres(nshg,nsd*nsd) 47 48c 49 integer ilwork(nlwork), 50 & invflx(nshg), nodflx(nshg) 51c 52 character*60 fname1, fmt1, fmt2, fnamer, fname2 53 character*13 fieldbflux 54 character*19 fieldwss 55 integer irstin, isize, nitems, ndofwss 56 integer iarray(50) ! integers read from headers 57 58 real*8, allocatable, dimension(:,:,:,:) :: xKebe, xGoC 59 integer, allocatable, dimension(:,:) :: ien2 60 integer, allocatable, dimension(:) :: map 61 real*8, allocatable, dimension(:,:) :: xmu2 62c 63c.... calculate the flux nodes 64c 65 numflx = 0 66 invflx = 0 67 nodflx = 0 68 do iblk = 1, nelblb 69 iel = lcblkb(1,iblk) 70 lcsyst = lcblkb(3,iblk) 71 nenl = lcblkb(5,iblk) 72 nenbl = lcblkb(6,iblk) 73 ndofl = lcblkb(8,iblk) 74 nshl = lcblkb(9,iblk) 75 nshlb = lcblkb(10,iblk) 76 npro = lcblkb(1,iblk+1) - iel 77 call flxNode (mienb(iblk)%p, miBCB(iblk)%p, invflx) 78 enddo 79c do i = 1, nshg 80 do i = 1, numnp 81 if (invflx(i) .ne. 0) then 82 numflx = numflx + 1 83 nodflx(numflx) = i 84 endif 85 enddo 86c 87c.... --------------------> interior elements <-------------------- 88c 89c.... initialize the arrays 90c 91 flxres = zero 92 flxLHS = zero 93 flxnrm = zero 94 if (numflx .ne. 0) then !we have flux nodes 95 qres = zero 96c 97c.... loop over the element-blocks 98c 99 lhs = 0 100 101 ires=2 ! shield e3ql from an unmapped lmassinv 102 ierrcalcsave=ierrcalc 103 ierrcalc=0 104 do iblk = 1, nelblk 105c 106c.... set up the parameters 107c 108 iel = lcblk(1,iblk) 109 nenl = lcblk(5,iblk) ! no. of vertices per element 110 nshl = lcblk(10,iblk) 111 ndofl = lcblk(8,iblk) 112 lcsyst = lcblk(3,iblk) 113 npro = lcblk(1,iblk+1) - iel 114 ngauss = nint(lcsyst) 115 allocate ( ien2(npro,nshl) ) 116 allocate ( xmu2(npro,maxsh)) 117 allocate ( map(npro) ) 118c 119c.... get the elements touching the boundary 120c 121 call mapConn( mien(iblk)%p, ien2, invflx, 122 & map, nshl, npro, 123 & npro2, nshg ) 124 125 nprold = npro 126 npro = npro2 127 128 if (npro .ne. 0) then 129 130 call mapArray( mxmudmi(iblk)%p, xmu2, map, 131 & maxsh, nprold) 132c 133c.... allocate the element matrices (though they're not needed) 134c 135 allocate ( xKebe(npro,9,nshl,nshl) ) 136 allocate ( xGoC (npro,4,nshl,nshl) ) 137c 138c.... compute and assemble the residuals 139c 140 call AsIGMR (y, ac, 141 & x, xmu2(1:npro,:), 142 & shp(lcsyst,1:nshl,:), 143 & shgl(lcsyst,:,1:nshl,:), 144 & ien2(1:npro,:), 145 & flxres, qres, 146 & xKebe, xGoC, 147 & rtmp) 148c 149 deallocate ( xKebe ) 150 deallocate ( xGoC ) 151 endif 152 deallocate ( ien2 ) 153 deallocate ( xmu2 ) 154 deallocate ( map ) 155c 156 enddo 157 ierrcalc=ierrcalcsave 158c 159c.... --------------------> boundary elements <-------------------- 160c 161 do iblk = 1, nelblb 162c 163c.... set up the parameters 164c 165 iel = lcblkb(1,iblk) 166 lcsyst = lcblkb(3,iblk) 167 nenl = lcblkb(5,iblk) 168 nshl = lcblkb(9,iblk) 169 nenbl = lcblkb(6,iblk) 170 nshlb = lcblkb(10,iblk) 171 npro = lcblkb(1,iblk+1) - iel 172 173 if(lcsyst.eq.3) lcsyst=nenbl 174c 175 if(lcsyst.eq.3 .or. lcsyst.eq.4) then 176 ngaussb = nintb(lcsyst) 177 else 178 ngaussb = nintb(lcsyst) 179 endif 180c 181c.... allocate the element matrices (though they're not needed) 182c 183 allocate ( xKebe(npro,9,nshl,nshl) ) 184c.... compute and assemble the residuals 185c 186 call AsBFlx (u, y, 187 & ac, x, 188 & shpb(lcsyst,1:nshl,:), 189 & shglb(lcsyst,:,1:nshl,:), 190 & mienb(iblk)%p, 191 & miBCB(iblk)%p, mBCB(iblk)%p, 192 & invflx, flxres, 193 & flxLHS, flxnrm, 194 & xKebe ) 195c 196 deallocate ( xKebe ) 197c 198c.... end of boundary element loop 199c 200 enddo 201c.... Communication needed before we take care of periodicity and 202c division of RHS by LHS ??? 203c 204 if(numpe > 1) then 205 call commu (flxres, ilwork, nflow, 'in ') 206 call commu (flxLHS, ilwork, 1 , 'in ') 207 call commu (flxnrm, ilwork, nsd , 'in ') 208 endif 209c 210c take care of periodic boundary conditions 211c 212 do j= 1,nshg 213 if ((btest(iBC(j),10))) then 214 i = iper(j) 215 flxLHS(i,1) = flxLHS(i,1) + flxLHS(j,1) 216 flxres(i,:) = flxres(i,:) + flxres(j,:) 217 endif 218 enddo 219 220 do j= 1,nshg 221 if ((btest(iBC(j),10))) then 222 i = iper(j) 223 flxLHS(j,1) = flxLHS(i,1) 224 flxres(j,:) = flxres(i,:) 225 endif 226 enddo 227c 228c call bc3per(iBC, flxres, iper, ilwork, nflow) 229 230c 231c.... integrated fluxes (aerodynamic forces update) 232c 233 flxTot = zero 234 do n = 1, numflx 235 flxTot = flxTot + flxres(nodflx(n),:) 236 enddo 237 Force(1) = flxTot(1) 238 Force(2) = flxTot(2) 239 Force(3) = flxTot(3) 240 else 241c need to call commu for procs. with no flux nodes 242 if(numpe > 1) then 243 call commu (flxres, ilwork, nflow, 'in ') 244 call commu (flxLHS, ilwork, 1 , 'in ') 245 call commu (flxnrm, ilwork, nsd , 'in ') 246 endif 247 endif ! make sure the zero numflux processors still commu 248c 249c.... only need to commu if we are going to print surface flux (or compute avg) since 250c the force calculation just sums flxres (and each on processor node 251c has his "piece" of the sum already). 252c 253 ntoutv=ntout 254 if ((irs .ge. 1) .and. (mod(lstep, ntoutv) .eq. 0) .or. 255 & istep .eq. nstep(itseq) .or. ioybar == 1) then 256 257 258c 259c need to zero the slaves to prevent counting twice 260c (actually unnecessary since flxres of boundary nodes will be counted n 261c times while flxlhs will be counted n times-> the ratio is still 262c correct 263c 264 265 ndofwss = 3; 266 wallssVec=rtmp(:,1:ndofwss) 267 268 if (numflx .eq. 0) then !no flux nodes 269 rtmp=zero 270 wallssVec = zero 271 272c need to call commu for procs. with no flux nodes 273 if(numpe > 1) then 274 call commu (flxres, ilwork, nflow, 'out') 275 call commu (flxLHS, ilwork, 1 , 'out') 276 call commu (flxnrm, ilwork, nsd , 'out') 277 endif 278 else 279c 280c.... ----------------------------> Solve <--------------------------- 281c 282c.... compute the viscous and heat fluxes 283c 284c 285c.... ----------------------------> Print <--------------------------- 286c 287c.... nodal fluxes 288c 289 do i = 1, 3 290 where ( (invflx .ne. 0) .and. (flxLHS(:,1) .ne. zero) ) 291 & flxres(:,i) = flxres(:,i) / flxLHS(:,1) 292 enddo 293c 294c.... normalize the outward normal 295c 296 temp = sqrt(flxnrm(:,1)**2 + flxnrm(:,2)**2 + flxnrm(:,3)**2) 297 where ( (invflx .ne. 0) .and. (temp .ne. zero) ) 298 flxnrm(:,1) = flxnrm(:,1) / temp 299 flxnrm(:,2) = flxnrm(:,2) / temp 300 flxnrm(:,3) = flxnrm(:,3) / temp 301 endwhere 302 303c 304c.... ----------------------------> Communications <------------------- 305c 306 if(numpe > 1) then 307 call commu (flxres, ilwork, nflow, 'out') 308 call commu (flxLHS, ilwork, 1 , 'out') 309 call commu (flxnrm, ilwork, nsd , 'out') 310 endif 311c 312 313 rtmp = zero 314 wallssVec = zero 315 316 do i=1, numnp 317 if (invflx(i) .ne. 0) then 318 rtmp(i,2:4) = flxres(i,1:3) !viscous flux 319c calculate the WSS 320 tn= flxres(i,1) * flxnrm(i,1) 321 & + flxres(i,2) * flxnrm(i,2) 322 & + flxres(i,3) * flxnrm(i,3) 323 324 wallssVec(i,1) = flxres(i,1) - tn * flxnrm(i,1) 325 wallssVec(i,2) = flxres(i,2) - tn * flxnrm(i,2) 326 wallssVec(i,3) = flxres(i,3) - tn * flxnrm(i,3) 327 328 endif 329 enddo 330 endif 331cc itmp = 1 332cc if (lstep .gt. 0) itmp = int(log10(float(lstep)))+1 333cc write (fmt1,"('(''flux.'',i',i1,',1x)')") itmp 334cc write (fname1,fmt1) lstep 335 336cc fname1 = trim(fname1) // cname2(myrank+1) 337 338cc open (unit=iflux, file=fname1, status='unknown', 339cc & form='formatted',err=997) 340 341c write (iflux) machin, nshg, lstep 342c write (iflux) rtmp(:,1:6) 343c 344c.... output the results 345c 346cc do n = 1, numflx 347cc k = nodflx(n) 348cc write (iflux,2000) k,invflx(k),flxLHS(k,1),(x(k,i),i=1,3), 349cc & (flxnrm(k,i), i=1,3), 350cc & (flxres(k,i), i=1,3) 351cc enddo 352cc close (iflux) 353c... output the results in the new format in restart.step#.proc# file 354 355cc itmp = 1 356cc if (lstep .gt. 0) itmp = int(log10(float(lstep)))+1 357cc write (fmt2,"('(''restart.'',i',i1,',1x)')") itmp 358cc write (fname2,fmt2) lstep 359 360cc fname2 = trim(fname2) // cname2(myrank+1) 361c 362c.... open input files 363c 364cc call openfile( fname2, 'append?', irstin ) 365 366cc fieldbflux = 'boundary flux' 367cc isize = nshg*ndof 368cc nitems = 3; 369cc iarray(1) = nshg 370cc iarray(2) = ndof 371cc iarray(3) = lstep 372cc call writeheader(irstin, fieldbflux,iarray, nitems, isize, 373cc & 'double', iotype ) 374 375c fieldbflux = 'boundary flux' 376cc nitems = nshg*ndof 377cc call writedatablock(irstin, fieldbflux,rtmp, nitems, 378cc & 'double', iotype) 379 380cc call closefile( irstin, "append" ) 381c call Write_boundaryflux(myrank,lstep,nshg,ndof,rtmp(:,1:ndof)) 382 383c wallss vectors into the restart file(s) 384 385 ntoutv=ntout 386 if ((irs .ge. 1) .and. (mod(lstep, ntoutv) .eq. 0) .or. 387 & istep .eq. nstep(itseq) ) then 388 389 call write_field(myrank,'a','wss',3, 390 & wallssVec,'d',nshg,ndofwss,lstep) 391 endif 392 393 endif 394c 395 return 396c 397c.... file error handling 398c 399997 call error ('bflux ','opening ', iflux) 400c 401c$$$1000 format(' ',a80,/,1x,i10,1p,3e20.7) 402 2000 format(i6,i6,10(2x,E12.5e2)) 403c$$$2001 format(1p,1x,i6,3e15.7) 404c 405c.... end 406c 407 end 408 409 410 subroutine flxNode(ienb, iBCB, flg) 411c--------------------------------------------------------------------- 412c 413c This routine flags the flux nodes 414c 415c---------------------------------------------------------------------- 416 include "common.h" 417c 418 integer flg(nshg), iBCB(npro,ndiBCB), 419 & ienb(npro, nshl), lnode(27) 420 421c 422c.... compute the nodes which lie on the boundary (hierarchic) 423c 424 call getbnodes(lnode) 425 426 do i=1, npro 427 if (nsrflist(iBCB(i,2)).eq.1) then 428 do j=1, nshlb 429 nodelcl = lnode(j) 430 flg(abs(ienb(i,nodelcl)))=flg(abs(ienb(i,nodelcl)))+1 431 enddo 432 endif 433 enddo 434c 435 return 436 end 437 438 439 subroutine mapConn( ien, ien2, mask, 440 & map, nshl, npro, 441 & npro2, nshg ) 442c----------------------------------------------------------------------- 443c 444c Create a condensed connectivity array based on the nodes in 445c mask. 446c 447c----------------------------------------------------------------------- 448 449 integer ien(npro,nshl), ien2(npro,nshl), mask(nshg), 450 & map(npro) 451 452 integer nshl, nshg, npro, npro2, i, iel 453 454c 455c.... first build the map 456c 457 map = 0 458 do i = 1, nshl 459 do iel = 1, npro 460 map(iel) = map(iel) + mask( abs(ien(iel,i)) ) 461 enddo 462 enddo 463 464 npro2 = 0 465 do iel = 1, npro 466 if ( map(iel) .gt. 0 ) then 467 npro2 = npro2 + 1 468 map(iel) = npro2 469 else 470 map(iel) = npro 471 endif 472 enddo 473c 474c.... create the condensed connectivity array 475c 476 if ( npro2 .gt. 0 ) then 477 do i = 1, nshl 478 do iel = 1, npro 479 ien2(map(iel),i) = ien(iel,i) 480 enddo 481 enddo 482 endif 483 484 return 485 end 486 487 488 subroutine mapArray( x, x2, map, 489 & nshl, nprold) 490c----------------------------------------------------------------------- 491c 492c Maps array x into array x2 based on the given map 493c 494c----------------------------------------------------------------------- 495 real*8 x(nprold,nshl), x2(nprold,nshl) 496 497 integer map(nprold) 498 499 integer nprold, nshl, i 500 501c 502c.... map the array 503c 504 do i = 1, nshl 505 x2(map(:),i) = x(:,i) 506 enddo 507 508 return 509 end 510