1 subroutine Bflux (y, ac, x, 2 & shp, shgl, shpb, 3 & shglb, nodflx, ilwork) 4c 5c---------------------------------------------------------------------- 6c 7c This routine : 8c 1. computes the boundary fluxes 9c 2. prints the results in the file [FLUX.lstep] 10c 3. prints the primitives variables in the files [OUTPUT.lstep] 11c and [OUTCHM.lstep] (2-D computations only) 12c 4. calls the restar routine 13c 14c 15c output: 16c in file FLUX.lstep: 17c run title 18c step number, run time 19c node_number 20c x_1 ... x_nsd ! coordinates 21c normal_1 ... normal_nsd ! outward normal direction 22c tau_1n ... tau_nsd n ! boundary viscous flux 23c q_n ! boundary head flux 24c ro u_1 ... u_nsd t p s Mach ! primitive variables 25c y_N2 y_O2 y_NO y_N y_O ! gas composition 26c 27c in file OUTPUT.lstep: 28c run title 29c step number, run time, gamma, Rgas 30c density, velocity vector and temperature 31c 32c in file OUTCHM.lstep: 33c run title 34c step number, run time, gamma, Rgas 35c pressure, entropy, Mach, y_N2, y_O2, y_NO, y_N and y_O 36c 37c 38c Zdenek Johan, Summer 1991. 39c---------------------------------------------------------------------- 40c 41 use pointer_data 42c 43 include "common.h" 44c 45 dimension y(nshg,ndof), ac(nshg,ndof), 46 & x(numnp,nsd), 47 & shp(MAXTOP,maxsh,MAXQPT), 48 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 49 & shpb(MAXTOP,maxsh,MAXQPT), 50 & shglb(MAXTOP,nsd,maxsh,MAXQPT), 51 & nodflx(numflx) 52c 53 dimension invflx(nshg), flxres(nshg,nflow), 54 & flxLHS(nshg,1), flxnrm(nshg,nsd), 55 & temp(nshg), itemp(numflx), 56 & q(nshg,ndof), qq(nshg,8), 57 & ilwork(nlwork) 58c 59 character(len=20) fname1, fname2, 60 & fmt1, fmt2 61c 62c.... ************************>> Restart <<*************************** 63c 64c.... set up the timer 65c 66 call timer ('Output ') 67c 68c.... call restart option 69c 70 q(:,1:ndof)=y(:,1:ndof) 71 72 if (irs .ge. 1) call restar ('out ', q,ac) 73c 74c No Flux output for now KENJ 75c 76 return 77c 78c.... scale the primitive variables for output 79c 80 q(:,1) = ro * q(:,1) 81 q(:,2) = vel * q(:,2) 82 q(:,3) = vel * q(:,3) 83 if (nsd .eq. 3) 84 & q(:,4) = vel * q(:,4) 85 q(:,nflow) = temper * q(:,nflow) 86 qq(:,1) = press * qq(:,1) 87 qq(:,2) = entrop * qq(:,2) 88c 89c.... *************************>> Output <<*************************** 90c 91c.... output only for 2-D computations 92c 93 if ((irs .eq. 2) .and. (nsd .eq. 2)) then 94c 95c.... get the file names ('output.'lstep) and ('outchm.'lstep) 96c 97 itmp = 1 98 if (lstep .gt. 0) itmp = int(log10(float(lstep)))+1 99 write (fmt1,"('(''output.'',i',i1,',1x)')") itmp 100 write (fmt2,"('(''outchm.'',i',i1,',1x)')") itmp 101 write (fname1,fmt1) lstep 102 write (fname2,fmt2) lstep 103c 104c.... open file 105c 106 if (ioform .eq. 0) then 107 open (unit=iout, file=fname1, status='unknown', err=995) 108 if (ichem .eq. 1) 109 & open (unit=ichmou, file=fname2, status='unknown', err=996) 110 else 111 open (unit=iout, file=fname1, status='unknown', err=995, 112 & form='unformatted') 113 if (ichem .eq. 1) 114 & open (unit=ichmou, file=fname2, status='unknown', err=996, 115 & form='unformatted') 116 endif 117c 118c.... output header 119c 120 if (ioform .eq. 0) then 121 write (iout, 1000) title, lstep, time, gamma, Rgas 122 if (ichem .eq. 1) 123 & write (ichmou,1000) title, lstep, time 124 else 125 write (iout) title 126 write (iout) nshg, lstep, time, gamma, Rgas 127c 128 if (ichem .eq. 1) then 129 write (ichmou) title 130 write (ichmou) nshg, lstep, time 131 endif 132 endif 133c 134 if (ioform .eq. 0) then 135 do n = 1, nshg 136 write (iout, 2000) n, (q(n,i), i=1,ndof) 137 if (ichem .eq. 1) 138 & write (ichmou,2000) n, (qq(n,i), i=1,8) 139 enddo 140 else 141 write (iout) q 142 if (ichem .eq. 1) write (ichmou) qq 143 endif 144c 145c.... close the files 146c 147 close (iout) 148 if (ichem .eq. 1) close (ichmou) 149c 150 endif 151c 152c.... stop the timer 153c 154 call timer ('Back ') 155c 156c.... *********************>> Boundary Fluxes <<********************** 157c 158 if ((irs .eq. 2) .and. (numflx .ne. 0)) then 159c 160c.... set up the timer 161c 162 call timer ('Bnd_Flux') 163c 164c.... calculate the inverse of nodflx 165c 166 itemp = nodflx 167c 168 invflx = 0 169 invflx(itemp) = (/ (i, i=1,numflx) /) 170c 171c.... --------------------> interior elements <-------------------- 172c 173c.... set up parameters 174c 175 intrul = intg (1,itseq) 176 intind = intpt (intrul) 177c 178 jump = min(iALE + iter-1, 1) 179 ires = 1 180 iprec = 0 181c 182c.... initialize the arrays 183c 184 flxres = zero 185 flxLHS = zero 186 flxnrm = zero 187c 188c.... loop over the element-blocks 189c 190 do iblk = 1, nelblk 191c 192c.... set up the parameters 193c 194c$$$ iel = lcblk(1,iblk) 195c$$$ nenl = lcblk(5,iblk) 196c$$$ mattyp = lcblk(7,iblk) 197c$$$ npro = lcblk(1,iblk+1) - iel 198c 199 nenl = lcblk(5,iblk) ! no. of vertices per element 200 iel = lcblk(1,iblk) 201 lelCat = lcblk(2,iblk) 202 lcsyst = lcblk(3,iblk) 203 iorder = lcblk(4,iblk) 204 nenl = lcblk(5,iblk) ! no. of vertices per element 205 nshl = lcblk(10,iblk) 206 mattyp = lcblk(7,iblk) 207 ndofl = lcblk(8,iblk) 208 nsymdl = lcblk(9,iblk) 209 npro = lcblk(1,iblk+1) - iel 210 ngauss = nint(lcsyst) 211c 212c 213 if (mattyp .ne. 0) cycle ! fluid elements only 214c 215c.... compute and assemble the residuals 216c 217 call AsIFlx (y, ac, 218 & x, mxmudmi(iblk)%p, 219 & shp, shgl, mien(iblk)%p, 220 & mmat(iblk)%p, flxres) 221c 222c.... end of interior element loop 223c 224 enddo 225c 226c.... --------------------> boundary elements <-------------------- 227c 228c.... loop over the elements 229c 230 do iblk = 1, nelblb 231c 232c.... set up the parameters 233c 234c$$$ iel = lcblkb(1,iblk) 235c$$$ nenl = 4 ! tetrahedral elements (4-vertices) 236c$$$ nenbl = 3 ! triangular boundary element faces 237c$$$ mattyp = lcblkb(7,iblk) 238c$$$ npro = lcblkb(1,iblk+1) - iel 239c 240c 241 iel = lcblkb(1,iblk) 242 lelCat = lcblkb(2,iblk) 243 lcsyst = lcblkb(3,iblk) 244 iorder = lcblkb(4,iblk) 245 nenl = lcblkb(5,iblk) ! no. of vertices per element 246 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 247 mattyp = lcblkb(7,iblk) 248 ndofl = lcblkb(8,iblk) 249 nshl = lcblkb(9,iblk) 250 nshlb = lcblkb(10,iblk) 251 npro = lcblkb(1,iblk+1) - iel 252 if(lcsyst.eq.3) lcsyst=nenbl 253 ngaussb = nintb(lcsyst) 254 255 if (mattyp .ne. 0) cycle ! fluid elements only 256c 257c.... compute and assemble the residuals 258c 259 call AsBFlx (y, x, 260 & shpb(lcsyst,1:nshl,:), 261 & shglb(lcsyst,:,1:nshl,:), 262 & mienb(iblk)%p, mmatb(iblk)%p, 263 & miBCB(iblk)%p, mBCB(iblk)%p, 264 & invflx, flxres, 265 & flxLHS, flxnrm) 266c 267c.... end of boundary element loop 268c 269 enddo 270c 271c.... ----------------------------> Communications <------------------- 272c 273 if(numpe > 1) then 274 call commu (flxres, ilwork, nflow, 'in ') 275 call commu (flxLHS, ilwork, 1 , 'in ') 276 call commu (flxnrm, ilwork, nsd , 'in ') 277 endif 278c 279c.... ----------------------------> Solve <--------------------------- 280c 281c.... compute the viscous and heat fluxes 282c 283 do i = 2, nflow 284 where ( (invflx .ne. 0) .and. (flxLHS(:,1) .ne. zero) ) 285 & flxres(:,i) = flxres(:,i) / flxLHS(:,1) 286 enddo 287c 288c.... normalize the outward normal 289c 290c if (nsd .eq. 2) then 291c 292c temp = sqrt(flxnrm(:,1)**2 + flxnrm(:,2)**2) 293c where ( (invflx .ne. 0) .and. (temp .ne. zero) ) 294c flxnrm(:,1) = flxnrm(:,1) / temp 295c flxnrm(:,2) = flxnrm(:,2) / temp 296c endwhere 297c 298c else 299c 300 temp = sqrt(flxnrm(:,1)**2 + flxnrm(:,2)**2 + flxnrm(:,3)**2) 301 where ( (invflx .ne. 0) .and. (temp .ne. zero) ) 302 flxnrm(:,1) = flxnrm(:,1) / temp 303 flxnrm(:,2) = flxnrm(:,2) / temp 304 flxnrm(:,3) = flxnrm(:,3) / temp 305 endwhere 306c 307c endif 308c 309c.... ----------------------------> Print <--------------------------- 310c 311c.... get the file name ('flux.'lstep) 312c 313 itmp = 1 314 if (lstep .gt. 0) itmp = int(log10(float(lstep)))+1 315 write (fmt1,"('(''flux.'',i',i1,',1x)')") itmp 316 write (fname1,fmt1) lstep 317c 318c.... open file 319c 320 open (unit=iflux, file=fname1, status='unknown', err=997) 321c 322c.... output the header 323c 324 write (iflux, 1000) title, lstep, time 325c 326c.... output the results 327c 328 do n = 1, numflx 329c 330 k = nodflx(n) 331 m = k ! global node number 332c 333 write (iflux, 2000) k, (x(m,i), i=1,nsd), 334 & (flxnrm(m,i), i=1,nsd), 335 & (flxres(m,i), i=2,nflow), 336 & (q(k,i), i=1,ndof), 337 & (qq(k,i), i=1,8) 338c 339 enddo 340c 341c.... close file 342c 343 close (iflux) 344c 345c.... stop the timer 346c 347 call timer ('Back ') 348c 349 endif 350c 351 return 352c 353c.... file error handling 354c 355995 call error ('bflux ','opening ', iout) 356996 call error ('bflux ','opening ', ichmou) 357997 call error ('bflux ','opening ', iflux) 358c 3591000 format(' ',a80,/,1x,i10,1p,3e20.7) 3602000 format(1p,1x,i6,23e15.7) 361c 362c.... end 363c 364 end 365