1c----------------------------------------------------------------------- 2c 3c module for time averaged statistics (conservative projection). 4c 5c----------------------------------------------------------------------- 6 module stats 7 8 integer nResDims, nSolDims, nLhsDims, nTimeStep, stsResFlg 9 integer stsCompFreq, stsWriteFreq, stsResetFreq, step1, 10 & stsType 11 12 real*8, allocatable :: stsVec(:,:) 13 14 real*8, allocatable :: stsReg(:) 15 real*8, allocatable :: stsMInv(:,:) 16 real*8, allocatable :: stsB(:,:) 17 real*8, allocatable :: stsDInv(:,:) 18 real*8, allocatable :: stsCInv(:,:) 19 20 real*8, allocatable :: stsPres(:), stsPresSqr(:), stsVel(:,:), 21 & stsVelSqr(:,:), stsVelReg(:,:), 22 & stsStress(:,:) 23 24 end module 25 26c----------------------------------------------------------------------- 27c create the new statistics arrays 28c----------------------------------------------------------------------- 29 subroutine initStats(x, iBC, iper, ilwork) 30 31 use stats 32 include "common.h" 33 34 real*8 x(numnp,3) 35 integer ilwork(nlwork), iper(nshg), iBC(nshg) 36 37 if (ipord .eq. 1) then 38 stsType = 1 39 else 40 stsType = 2 41 endif 42 43 stsWriteFreq = 200 44 45 nResDims = 11 46 nSolDims = 10 47 nLhsDims = 19 48 49 allocate ( stsVec(nshg,nResDims) ) 50 51 if (stsType .eq. 1) then 52 allocate ( stsReg(nshg) ) 53 allocate ( stsMInv(nshg,6) ) 54 allocate ( stsB(nshg,3) ) 55 allocate ( stsDInv(nshg,3) ) 56 allocate ( stsCInv(nshg,6) ) 57 endif 58 59 allocate ( stsPres(nshg) ) 60 allocate ( stsPresSqr(nshg) ) 61 allocate ( stsVel(nshg,3) ) 62 allocate ( stsVelSqr(nshg,6) ) 63 allocate ( stsVelReg(nshg,3) ) 64 allocate ( stsStress(nshg,6) ) 65 66 stsPres = 0.0 67 stsPresSqr = 0.0 68 stsVel = 0.0 69 stsVelSqr = 0.0 70 stsVelReg = 0.0 71 stsStress = 0.0 72 73 step1 = lstep+1 74 nTimeStep = 0 75 stsResFlg = 0 76 77 if (stsType .eq. 1) then 78 call elmStatsLhs( x, iBC, iper, ilwork ) 79 call stsInitLhs( nshg ) 80 endif 81 82 return 83 end 84 85c----------------------------------------------------------------------- 86c compute the Lhs matrices needed for the conservative projection 87c of the statistics 88c----------------------------------------------------------------------- 89 subroutine stsInitLhs(nshg) 90 91 use stats 92 integer nshg 93 94 real*8 res(nResDims), reg, det, r0, r1, r2, r3, r4, r5, 95 & d0, d1, d2, c0, c1, c2, c3, c4, c5 96 integer i 97 98c 99c.... build the regularization 100c 101 do i = 1, nshg 102 res = stsVec(i,:) 103 reg = res(1) * res(2) * res(3) 104 det = res(1) * (res(2) * res(3) - res(5) * res(5)) 105 & + res(4) * (res(5) * res(6) - res(4) * res(3)) 106 & + res(6) * (res(4) * res(5) - res(2) * res(6)) 107 108 if ( det .gt. 1.d-10*reg .and. reg .ne. 0 ) then 109 stsReg(i) = 0 110 else 111 stsReg(i) = (res(1) + res(2) + res(3)) / 1000. 112 endif 113 enddo 114 115c 116c.... form M and factorize 117c 118 do i = 1, nshg 119 res = stsVec(i,:) 120 reg = stsReg(i) 121 r0 = res(1) + reg 122 r1 = res(2) + reg 123 r2 = res(3) + reg 124 r3 = res(4) 125 r4 = res(5) 126 r5 = res(6) 127 128 det = r0 * (r1 * r2 - r4 * r4) 129 & + r3 * (r4 * r5 - r3 * r2) 130 & + r5 * (r3 * r4 - r1 * r5) 131 det = 1.0/det 132 133 stsMInv(i,1) = det * (r1 * r2 - r4 * r4) 134 stsMInv(i,2) = det * (r0 * r2 - r5 * r5) 135 stsMInv(i,3) = det * (r0 * r1 - r3 * r3) 136 stsMInv(i,4) = det * (r4 * r5 - r2 * r3) 137 stsMInv(i,5) = det * (r3 * r5 - r0 * r4) 138 stsMInv(i,6) = det * (r3 * r4 - r1 * r5) 139 enddo 140 141c 142c.... form B, DInv and CInv 143c 144 do i = 1, nshg 145 res = stsVec(i,:) 146 reg = stsReg(i) 147 r0 = res(1) 148 r1 = res(2) 149 r2 = res(3) 150 r3 = res(4) 151 r4 = res(5) 152 r5 = res(6) 153 d0 = 1. / ( reg + r0 ) 154 d1 = 1. / ( reg + r1 ) 155 d2 = 1. / ( reg + r2 ) 156 stsDInv(i,1) = d0 157 stsDInv(i,2) = d1 158 stsDInv(i,3) = d2 159 stsB(i,1) = r3 160 stsB(i,2) = r4 161 stsB(i,3) = r5 162 c0 = r0 + r1 - r3 * r3 * (d0 + d1) + reg 163 c1 = r1 + r2 - r4 * r4 * (d1 + d2) + reg 164 c2 = r2 + r0 - r5 * r5 * (d2 + d0) + reg 165 c3 = r5 - r3 * r4 * d1 166 c4 = r3 - r4 * r5 * d2 167 c5 = r4 - r5 * r3 * d0 168 det = c0 * (c1 * c2 - c4 * c4) 169 & + c3 * (c4 * c5 - c3 * c2) 170 & + c5 * (c3 * c4 - c1 * c5) 171 det = 1. / det 172 stsCInv(i,1) = det * (c1 * c2 - c4 * c4) 173 stsCInv(i,2) = det * (c0 * c2 - c5 * c5) 174 stsCInv(i,3) = det * (c0 * c1 - c3 * c3) 175 stsCInv(i,4) = det * (c4 * c5 - c2 * c3) 176 stsCInv(i,5) = det * (c3 * c5 - c0 * c4) 177 stsCInv(i,6) = det * (c3 * c4 - c1 * c5) 178 enddo 179 180 return 181 end 182 183 184c----------------------------------------------------------------------- 185c collect the desired statistics 186c----------------------------------------------------------------------- 187 subroutine stsGetStats( y, yold, ac, acold, 188 & u, uold, x, 189 & shp, shgl, shpb, shglb, 190 & iBC, BC, iper, ilwork, 191 & rowp, colm, lhsK, lhsP ) 192 193 use stats 194 195 include "common.h" 196 197 real*8 y(nshg,ndof), yold(nshg,ndof), 198 & ac(nshg,ndof), acold(nshg,ndof), 199 & u(nshg,nsd), uold(nshg,nsd), 200 & shp(MAXTOP,maxsh,MAXQPT), shgl(MAXTOP,nsd,maxsh,MAXQPT), 201 & shpb(MAXTOP,maxsh,MAXQPT), 202 & shglb(MAXTOP,nsd,maxsh,MAXQPT), 203 & BC(nshg,ndofBC), lhsK(9,nnz_tot), 204 & lhsP(4,nnz_tot), x(numnp,nsd) 205 206 integer iBC(nshg), iper(nshg), 207 & ilwork(nlwork), rowp(nshg*nnz), 208 & colm(nshg+1) 209 210 211 real*8 yAlpha(nshg,ndof), acAlpha(nshg,ndof), 212 & uAlpha(nshg,nsd), 213 & res(nResDims), MInv(6), 214 & DInv(3), B(3), 215 & CInv(6) 216 217 real*8 u1, u2, u3, r0, r1, r2, r3, r4, r5, t3, t4, t5 218 219 integer i 220 221 nTimeStep = nTimeStep + 1 222c 223c.... compute solution at n+alpha 224c 225 call itrYAlpha( uold, yold, acold, 226 & u, y, ac, 227 & uAlpha, yAlpha, acAlpha) 228 229c 230c.... assemble the residual 231c 232 if (stsType .eq. 1) then 233 call elmStatsRes( yAlpha, acAlpha, x, shp, shgl, 234 & shpb, shglb, iBC, BC, 235 & iper, ilwork, rowp, colm, 236 & lhsK, lhsP ) 237 238c 239c.... compute the statistics 240c 241 do i = 1, nshg 242 res = stsVec(i,:) 243 reg = stsReg(i) 244 245 MInv = stsMInv(i,:) 246 DInv = stsDInv(i,:) 247 B = stsB(i,:) 248 CInv = stsCInv(i,:) 249 250 u1 = yAlpha(i,1) 251 u2 = yAlpha(i,2) 252 u3 = yAlpha(i,3) 253 254 stsPres(i) = stsPres(i) + y(i,4) 255 stsPresSqr(i) = stsPresSqr(i) + y(i,4)*y(i,4) 256 257 r0 = res(1) + reg * u1 258 r1 = res(2) + reg * u2 259 r2 = res(3) + reg * u3 260 261 stsVel(i,1) = stsVel(i,1) 262 & + MInv(1) * r0 + MInv(4) * r1 + MInv(6) * r2 263 stsVel(i,2) = stsVel(i,2) 264 & + MInv(4) * r0 + MInv(2) * r1 + MInv(5) * r2 265 stsVel(i,3) = stsVel(i,3) 266 & + MInv(6) * r0 + MInv(5) * r1 + MInv(3) * r2 267 268 stsVelReg(i,1) = stsVelReg(i,1) + u1 269 stsVelReg(i,2) = stsVelReg(i,2) + u2 270 stsVelReg(i,3) = stsVelReg(i,3) + u3 271 272 r0 = res(1) * u1 + reg * u1 * u1 273 r1 = res(2) * u2 + reg * u2 * u2 274 r2 = res(3) * u3 + reg * u3 * u3 275 r3 = res(1) * u2 + res(2) * u1 + reg * u2 * u1 276 r4 = res(2) * u3 + res(3) * u2 + reg * u3 * u2 277 r5 = res(3) * u1 + res(1) * u3 + reg * u1 * u3 278 r0 = DInv(1) * r0 279 r1 = DInv(2) * r1 280 r2 = DInv(3) * r2 281 r3 = r3 - B(1) * (r0 + r1) 282 r4 = r4 - B(2) * (r1 + r2) 283 r5 = r5 - B(3) * (r2 + r0) 284 t3 = CInv(1) * r3 + CInv(4) * r4 + CInv(6) * r5 285 t4 = CInv(4) * r3 + CInv(2) * r4 + CInv(5) * r5 286 t5 = CInv(6) * r3 + CInv(5) * r4 + CInv(3) * r5 287 288 stsVelSqr(i,1) = stsVelSqr(i,1) 289 & + r0 - DInv(1) * (B(1) * t3 + B(3) * t5) 290 stsVelSqr(i,2) = stsVelSqr(i,2) 291 & + r1 - DInv(2) * (B(2) * t4 + B(1) * t3) 292 stsVelSqr(i,3) = stsVelSqr(i,3) 293 & + r2 - DInv(3) * (B(3) * t5 + B(2) * t4) 294 295 stsVelSqr(i,4) = stsVelSqr(i,4) + t3 296 stsVelSqr(i,5) = stsVelSqr(i,5) + t4 297 stsVelSqr(i,6) = stsVelSqr(i,6) + t5 298 299 r0 = res(4) 300 r1 = res(5) 301 r2 = res(6) 302 r3 = res(7) 303 r4 = res(8) 304 r5 = res(9) 305 306 r0 = DInv(1) * r0 307 r1 = DInv(2) * r1 308 r2 = DInv(3) * r2 309 310 r3 = r3 - B(1) * (r0 + r1) 311 r4 = r4 - B(2) * (r1 + r2) 312 r5 = r5 - B(3) * (r2 + r0) 313 314 t3 = CInv(1) * r3 + CInv(4) * r4 + CInv(6) * r5 315 t4 = CInv(4) * r3 + CInv(2) * r4 + CInv(5) * r5 316 t5 = CInv(6) * r3 + CInv(5) * r4 + CInv(3) * r5 317 318 stsStress(i,1) = stsStress(i,1) 319 & + r0 - DInv(1) * (B(1) * t3 + B(3) * t5) 320 stsStress(i,2) = stsStress(i,2) 321 & + r1 - DInv(2) * (B(2) * t4 + B(1) * t3) 322 stsStress(i,3) = stsStress(i,3) 323 & + r2 - DInv(3) * (B(3) * t5 + B(2) * t4) 324 stsStress(i,4) = stsStress(i,4) + t3 325 stsStress(i,5) = stsStress(i,5) + t4 326 stsStress(i,6) = stsStress(i,6) + t5 327 enddo 328 else if (stsType .eq. 2) then 329 330 call evalAtInterp( yAlpha, stsVec, x, 331 & nResDims, nshape) 332 333 do i = 1, nshg 334 335 u1 = stsVec(i,1) 336 u2 = stsVec(i,2) 337 u3 = stsVec(i,3) 338 339 stsPres(i) = stsPres(i) + stsVec(i,4) 340 stsPresSqr(i) = stsPresSqr(i) + stsVec(i,4)*stsVec(i,4) 341 342 stsVel(i,1) = stsVel(i,1) + u1 343 stsVel(i,2) = stsVel(i,2) + u2 344 stsVel(i,3) = stsVel(i,3) + u3 345 346 stsVelSqr(i,1) = stsVelSqr(i,1) + u1*u1 347 stsVelSqr(i,2) = stsVelSqr(i,2) + u2*u2 348 stsVelSqr(i,3) = stsVelSqr(i,3) + u3*u3 349 stsVelSqr(i,4) = stsVelSqr(i,4) + u1*u2 350 stsVelSqr(i,5) = stsVelSqr(i,5) + u2*u3 351 stsVelSqr(i,6) = stsVelSqr(i,6) + u3*u1 352 353 stsStress(i,1) = stsStress(i,1) + stsVec(i,6) 354 stsStress(i,2) = stsStress(i,2) + stsVec(i,7) 355 stsStress(i,3) = stsStress(i,3) + stsVec(i,8) 356 stsStress(i,4) = stsStress(i,4) + stsVec(i,9) 357 stsStress(i,5) = stsStress(i,5) + stsVec(i,10) 358 stsStress(i,6) = stsStress(i,6) + stsVec(i,11) 359 360 enddo 361 endif 362 363 if ( mod(nTimeStep,stsWriteFreq) .eq. 0 .or. 364 & nTimeStep .eq. nstep(itseq) ) then 365 call stsWriteStats() 366 endif 367 368c$$$ if ( mod( nTimeStep, stsResetFreq) .eq. 0 ) then 369c$$$ stsPres = 0.0 370c$$$ stsPresSqr = 0.0 371c$$$ stsVel = 0.0 372c$$$ stsVelSqr = 0.0 373c$$$ stsVelReg = 0.0 374c$$$ stsStress = 0.0 375c$$$ 376c$$$ nTimeStep = 0 377c$$$ endif 378 379 return 380 end 381 382c----------------------------------------------------------------------- 383c collect the desired statistics 384c----------------------------------------------------------------------- 385 subroutine stsWriteStats() 386 387 use stats 388 include "common.h" 389 390 integer iofile, step2, itmp1, itmp,iofile2 391 character*30 fname, fmt1 392 character*5 cname 393 character*1 dash 394 real*8 outvec(nshg,19) 395c 396c.... open the output file 397c 398 iofile = 39 399 step2 = lstep+1 ! current time step 400 itmp = 1 401 itmp1 = 1 402 if (step1 .gt. 0) itmp = int(log10(float(step1)))+1 403 if (step2 .gt. 0) itmp1 = int(log10(float(step2)))+1 404 dash = '-' 405 write (fmt1, 406 & "('(''stats.'',i',i1,',a1,i',i1,',1x)')") 407 & itmp,itmp1 408 write (fname,fmt1) step1,dash,step2 409 fname = trim(fname) // cname(myrank+1) 410 open ( unit = iofile, file = fname, status = 'unknown', 411 & form = 'unformatted') 412c 413c.... write the statistics 414c 415 outvec(:,1) = stsPres(:) 416 outvec(:,2:4) = stsVel(:,:) 417c outvec(:,2:4) = stsVelReg(:,:) 418 outvec(:,5) = zero ! later will be temperature 419 outvec(:,6) = stsPresSqr(:) 420 outvec(:,7:12) = stsVelSqr(:,:) 421 outvec(:,13) = zero ! later wil be tempSqr 422 outvec(:,14:19) = stsStress(:,:) 423 424 write (iofile) numnp, nshg, nTimeStep 425 write (iofile) outvec(1:nshg,:) 426 close (iofile) 427 428c$$$ write (iofile) numnp, numnp, nTimeStep 429c$$$ write (iofile) outvec(1:numnp,:) 430c$$$ close (iofile) 431 432 iofile2 = 40 433c$$$ open (unit=iofile2,file='stats.asc',status='unknown') 434c$$$c$$$c 435c$$$c$$$c.... pressure, velocity, temperature 436c$$$c$$$c 437c$$$c$$$ write (iofile2) outvec(:,1:5) 438c$$$c$$$c 439c$$$c$$$c.... pressSqr, u_i u_j, tempSqr 440c$$$c$$$c 441c$$$c$$$ write (iofile2) outvec(:,6:13) 442c$$$c$$$c 443c$$$c$$$c.... viscous stress 444c$$$c$$$c 445c$$$c$$$ write (iofile2) outvec(:,14:19) 446c$$$c$$$ 447c$$$ 448c$$$c 449c$$$c.... write the statistics 450c$$$c 451c$$$ write (iofile2,*) 'nNodes ',numnp 452c$$$ write (iofile2,*) 'power ',2 453c$$$ write (iofile2,*) 'nTimeSteps = ', nTimeStep 454c$$$c 455c$$$c.... velocity 456c$$$c 457c$$$ write (iofile2,*) 'vel 3 ',numnp 458c$$$ do i = 1, numnp 459c$$$ write (iofile2,111) stsVel(i,1), stsVel(i,2), stsVel(i,3) 460c$$$ enddo 461c$$$c 462c$$$c.... pressure 463c$$$c 464c$$$ write (iofile2,*) 'pres 1 ',numnp 465c$$$ do i = 1, numnp 466c$$$ write (iofile2,112) stsPres(i) 467c$$$ enddo 468c$$$c 469c$$$c.... velSqr 470c$$$c 471c$$$ write (iofile2,*) 'velSqr 6 ',numnp 472c$$$ do i = 1, numnp 473c$$$ write (iofile2,113) (stsVelSqr(i,j),j=1,6) 474c$$$ enddo 475c$$$c 476c$$$c.... presSqr 477c$$$c 478c$$$ write (iofile2,*) 'presSqr 1 ',numnp 479c$$$ do i = 1, numnp 480c$$$ write (iofile2,112) stsPresSqr(i) 481c$$$ enddo 482c$$$c 483c$$$c.... stress 484c$$$c 485c$$$ write (iofile2,*) 'stress 6 ',numnp 486c$$$ do i = 1, numnp 487c$$$ write (iofile2,113) (stsStress(i,j),j=1,6) 488c$$$ enddo 489c$$$c 490c$$$c.... velocity 491c$$$c 492c$$$ write (iofile2,*) 'vel 3 ',numnp 493c$$$ do i = 1, numnp 494c$$$ write (iofile2,111) (stsVelReg(i,j),j=1,3) 495c$$$ enddo 496c$$$ 497c$$$ close (iofile2) 498 499 111 format(1p,3e24.16) 500 112 format(1p, e24.16) 501 113 format(1p,6e24.16) 502 503 return 504 end 505