1 module aveall 2 3 real*8, allocatable :: ylist(:) 4 integer, allocatable :: ifathe(:,:) 5 6 end module 7 8c------------------------------------------------------------------------------ 9 10 subroutine setave 11 12 use aveall 13 14 include "common.h" 15 16 idim = numel*intg(1,1) 17 allocate ( ylist(idim) ) 18 allocate ( ifathe(numel,maxsh) ) 19 20c.... return 21 22 return 23 end 24 25c------------------------------------------------------------------------------ 26 27 subroutine aveprep(shp,x) 28 29 use pointer_data 30 use aveall 31 32 include "common.h" 33 34 dimension shp(MAXTOP,maxsh,MAXQPT) 35 dimension x(numnp,3) 36 integer nlist 37 38 nlist = 0 39 ylist = zero 40 lfathe = 0 41 do iblk = 1,nelblk 42 lcsyst = lcblk(3,iblk) 43 iel = lcblk(1,iblk) !Element number where this block begins 44 npro = lcblk(1,iblk+1) - iel 45 lelCat = lcblk(2,iblk) 46 nenl = lcblk(5,iblk) 47 nshl = lcblk(10,iblk) 48 inum = iel + npro - 1 49 50 ngauss = nint(lcsyst) 51 52 call getylist( ylist, ifathe(iel:inum,:), shp(lcsyst,1:nshl,:), 53 & x, mien(iblk)%p) 54 55 56 enddo 57 58 return 59 end 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 subroutine getnu (ien, strl, xmudmi, cdelsq, lfathe) 76 77 include "common.h" 78 79 dimension ien(npro,nshl), strl(npro,ngauss), 80 & xmudmi(npro,ngauss), cdelsq(nlist), 81 & lfathe(npro,ngauss) 82 83 do intp = 1, ngauss 84 do iel = 1, npro 85 86 xmudmi(iel,intp) = cdelsq( lfathe(iel,intp) )* 87 & strl(iel,intp) 88 89 enddo 90 enddo 91 92c 93c local clipping 94c 95 rmu=datmat(1,2,1) 96cc xmudmi=min(xmudmi,1000.0*rmu) !don't let it get larger than 1000 mu 97cc xmudmi=max(xmudmi, -rmu) ! don't let (xmudmi + mu) < 0 98 99 do int = 1, ngauss 100 xmudmi(:,int) = 101 & min(xmudmi(:,int),1000.0*rmu) !don't let it get larger than 1000 mu 102 xmudmi(:,int) = 103 & max(xmudmi(:,int), -rmu) ! don't let (xmudmi + mu) < 0 104 enddo 105 106 return 107 end 108 subroutine getxnut (fres, xden, xnum, ien, shp) 109 110 include "common.h" 111 112 dimension ien(npro,nshl), fres(nshg,22), 113 & shp(nshl,maxsh), fresli(npro,22), 114 & fresl(npro,nshl,22), strnrmi(npro), 115 & xnutl(npro), fwr(npro), 116 & xmij(npro,6), xlij(npro,6), 117 & xdenli(npro), xnumli(npro), 118 & xdenl(npro), xnuml(npro), 119 & xden(1), xnum(1) 120 121 122 call local (fres, fresl, ien, 22, 'gather ') 123 124 xnuml = zero 125 xdenl = zero 126 127 do intp = 1, ngauss 128 129 fresli = zero 130 do i = 1, nenl 131 do j = 1, 22 132 fresli(:,j) = fresli(:,j) + shp(i,intp)*fresl(:,i,j) 133 enddo 134 enddo 135 136 strnrmi(:) = sqrt( 137 & two * (fresli(:,10)**2 + fresli(:,11)**2 + fresli(:,12)**2) 138 & + four * ( fresli(:,13)**2 + fresli(:,14)**2 + 139 & fresli(:,15)**2 ) ) 140 141 fwr = fwr1 * fresli(:,22) * strnrmi(:) 142 143 xmij(:,1) = -fwr 144 & * fresli(:,10) + fresli(:,16) 145 146 xmij(:,2) = -fwr 147 & * fresli(:,11) + fresli(:,17) 148 149 xmij(:,3) = -fwr 150 & * fresli(:,12) + fresli(:,18) 151 152 xmij(:,4) = -fwr * fresli(:,13) + fresli(:,19) 153 xmij(:,5) = -fwr * fresli(:,14) + fresli(:,20) 154 xmij(:,6) = -fwr * fresli(:,15) + fresli(:,21) 155 156 fresli(:,22) = one/fresli(:,22) 157 158 xlij(:,1) = fresli(:,4) - fresli(:,1) * 159 & fresli(:,1) * fresli(:,22) 160 xlij(:,2) = fresli(:,5) - fresli(:,2) * 161 & fresli(:,2) * fresli(:,22) 162 xlij(:,3) = fresli(:,6) - fresli(:,3) * 163 & fresli(:,3) * fresli(:,22) 164 xlij(:,4) = fresli(:,7) - fresli(:,1) * 165 & fresli(:,2) * fresli(:,22) 166 xlij(:,5) = fresli(:,8) - fresli(:,1) * 167 & fresli(:,3) * fresli(:,22) 168 xlij(:,6) = fresli(:,9) - fresli(:,2) * 169 & fresli(:,3) * fresli(:,22) 170 171 xnumli(:) = xlij(:,1) * xmij(:,1) 172 & + xlij(:,2) * xmij(:,2) 173 & + xlij(:,3) * xmij(:,3) 174 & + two * (xlij(:,4) * xmij(:,4) 175 & + xlij(:,5) * xmij(:,5) 176 & + xlij(:,6) * xmij(:,6)) 177 xdenli(:) = xmij(:,1) * xmij(:,1) 178 & + xmij(:,2) * xmij(:,2) 179 & + xmij(:,3) * xmij(:,3) 180 & + two * (xmij(:,4) * xmij(:,4) 181 & + xmij(:,5) * xmij(:,5) 182 & + xmij(:,6) * xmij(:,6)) 183 184c xdenli = xdenli 185 186 xdenl(:) = xdenl(:) + xdenli(:) 187 xnuml(:) = xnuml(:) + xnumli(:) 188 189 enddo 190 191 do nel = 1, npro 192 xden(1) = xden(1) + xdenl(nel) 193 xnum(1) = xnum(1) + xnuml(nel) 194 enddo 195 196 return 197 198 end 199 200 201 202 203 204 205 206 207 208 subroutine getylist (ylist, lfathe, shp, x, ien) 209 210 include "common.h" 211 212 dimension lfathe(npro,maxsh), shp(nshl,maxsh), 213 & x(numnp,3), ien(npro,nshl), 214 & xl(npro,nenl,nsd), ylist(idim) 215 216 integer intp 217 218 call localx (x, xl, ien, 3, 'gather ') 219 220 221 do nel = 1, npro 222 do intp = 1, ngauss 223 224c... Compute the y-coordinate of the current quad. pt. and call it yint. 225 226 yint = zero 227 xint = zero 228 zint = zero 229 do n = 1, nenl 230 yint = yint + xl(nel,n,2)*shp(n,intp) 231 xint = xint + xl(nel,n,1)*shp(n,intp) 232 zint = zint + xl(nel,n,3)*shp(n,intp) 233 enddo 234 235c... Check ylist to see if yint is already included in ylist. 236 237 if(nlist .eq. 0)then ! In this case yint is definitely not in ylist. 238 nlist = nlist + 1 239 ylist(nlist) = yint 240 lfathe(nel,intp) = nlist 241 242 else 243 244 do ilist = 1, nlist 245 imatch = ilist 246 if( abs(yint-ylist(ilist)) .lt. 0.00001 ) then 247 lfathe(nel,intp) = imatch 248 goto 5 249 endif 250 enddo 251 252 nlist = nlist + 1 253 ylist(nlist) = yint 254 lfathe(nel,intp) = nlist 255 256 endif 257 258 5 yint = zero 259 260c if( nlist .eq. 5)then 261c write(*,*)'ylist yint=',ylist(4),ylist(5) 262c stop 263c endif 264 265 enddo ! End loop over quad. pts. 266 enddo ! End loop over elements in current block. 267 268 return 269 end 270 subroutine projdmc (y, shgl, shp, 271 & iper, ilwork, x) 272 273 use pointer_data 274 275 use aveall ! This module helps us average cdelsq computed at quad pts 276 277 use quadfilt ! This module gives us shglf(maxtp,nsd,maxsh,ngaussf), 278c shpf(maxtp,maxsh,ngaussf), and Qwtf(maxtp,ngaussf). 279c Shpf and shglf are the shape funciotns and their 280c gradient evaluated using the quadrature rule desired 281c for computing the dmod. Qwt contains the weights of the 282c quad. points. 283 284 include "common.h" 285 include "mpif.h" 286 include "auxmpi.h" 287 288c 289 dimension fres(nshg,24), fwr(nshg), 290 & strnrm(nshg), cdelsq(nshg), 291 & strl(numel,ngauss), 292 & y(nshg,5), yold(nshg,5), 293 & iper(nshg), 294 & ilwork(nlwork), 295 & x(numnp,3), 296 & shgl(MAXTOP,nsd,maxsh,MAXQPT), shp(MAXTOP,maxsh,MAXQPT), 297 & xden(idim), xnum(idim), 298 & xdent(idim), xnumt(idim) 299 300 fres = zero 301 yold(:,1)=y(:,4) 302 yold(:,2:4)=y(:,1:3) 303c 304c hack in an interesting velocity field (uncomment to test dmod) 305c 306c yold(:,5) = 1.0 307c yold(:,2) = 2.0*x(:,1) - 3*x(:,2) 308c yold(:,3) = 3.0*x(:,1) + 4.0*x(:,2) 309c yold(:,4) = 4.0*x(:,1) + x(:,2) + x(:,3) 310c yold(:,1) = Rgas * yold(:,5) ! Necessary to make model suitable 311c suitable for the 312 313 314 intrul=intg(1,itseq) 315 intind=intpt(intrul) 316 317 do iblk = 1,nelblk 318 lcsyst = lcblk(3,iblk) 319 iel = lcblk(1,iblk) !Element number where this block begins 320 npro = lcblk(1,iblk+1) - iel 321 lelCat = lcblk(2,iblk) 322 nenl = lcblk(5,iblk) 323 nshl = lcblk(10,iblk) 324 inum = iel + npro - 1 325 326 ngauss = nint(lcsyst) 327 ngaussf = nintf(lcsyst) 328 329 call asithf (yold, x, strl(iel:inum,:), mien(iblk)%p, fres, 330 & shglf(lcsyst,:,1:nshl,:), 331 & shpf(lcsyst,1:nshl,:),Qwtf(lcsyst,1:ngaussf)) 332 333 enddo 334c 335 336 if(numpe>1) call commu (fres, ilwork, 24, 'in ') 337c 338c account for periodicity in filtered variables 339c 340 do j = 1,nshg 341 i = iper(j) 342 if (i .ne. j) then 343 fres(i,:) = fres(i,:) + fres(j,:) 344 endif 345 enddo 346 347 do j = 1,nshg 348 i = iper(j) 349 if (i .ne. j) then 350 fres(j,:) = fres(i,:) 351 endif 352 enddo 353 354 if(numpe>1)then 355 call commu (fres, ilwork, 24, 'out') 356 endif 357 358 fres(:,23) = one / fres(:,23) 359 do j = 1,22 360 fres(:,j) = fres(:,j) * fres(:,23) 361 enddo 362 363 xden = zero 364 xdent = zero 365 xnum = zero 366 xnumt = zero 367 368 do iblk = 1,nelblk 369 lcsyst = lcblk(3,iblk) 370 iel = lcblk(1,iblk) !Element number where this block begins 371 npro = lcblk(1,iblk+1) - iel 372 lelCat = lcblk(2,iblk) 373 nenl = lcblk(5,iblk) 374 inum = iel + npro - 1 375 ngauss = nint(lcsyst) 376 call smagcoeff (fres(:,1:22), xden, xnum, mien(iblk)%p, 377 & shp(lcsyst,1:nshl,:), ifathe(iel:inum,:)) 378 379 enddo 380 381 if(numpe.gt.1)then 382 call drvAllreduce( xden, xdent, idim ) 383 call drvAllreduce( xnum, xnumt, idim ) 384 else 385 xdent = xden 386 xnumt = xnum 387 endif 388 389 do i = 1, nlist 390 cdelsq(i) = pt5*xnumt(i)/xdent(i) 391 enddo 392 393c... Uncomment for averaging over all directions 394 395c sumc1 = 0.0 396c sumc2 = 0.0 397c do i = 1, nlist 398c sumc1 = sumc1 + xnumt(i) 399c sumc2 = sumc2 + xdent(i) 400c enddo 401c xfact = pt5*sumc1/sumc2 402c cdelsq(:) = xfact 403c if(myrank .eq. master)then 404c write(22,*)cdelsq(1) 405c endif 406 407c... End of averaging over all directions. 408 409 do iblk = 1,nelblk 410 lcsyst = lcblk(3,iblk) 411 iel = lcblk(1,iblk) !Element number where this block begins 412 npro = lcblk(1,iblk+1) - iel 413 lelCat = lcblk(2,iblk) 414 nenl = lcblk(5,iblk) 415 nshl = lcblk(10,iblk) 416 inum = iel + npro - 1 417 418 ngauss = nint(lcsyst) 419 ngaussf = nintf(lcsyst) 420 421 if (ngauss .ne. ngaussf) then 422 call getstrl (yold, x, mien(iblk)%p, 423 & strl(iel:inum,:), shgl(lcsyst,:,1:nshl,:), 424 & shp(lcsyst,1:nshl,:)) 425 endif 426 427 enddo 428 429 do iblk = 1,nelblk 430 lcsyst = lcblk(3,iblk) 431 iel = lcblk(1,iblk) 432 npro = lcblk(1,iblk+1) - iel 433 lelCat = lcblk(2,iblk) 434 inum = iel + npro - 1 435 436 ngauss = nint(lcsyst) 437 438 call getnu (mien(iblk)%p, strl(iel:inum,:), 439 & mxmudmi(iblk)%p,cdelsq,ifathe(iel:inum,:)) 440 enddo 441 442 443 return 444 end 445 446 447 448 subroutine smagcoeff (fres, xden, xnum, ien, shp, lfathe) 449 450 include "common.h" 451 452 dimension ien(npro,nshl), fres(nshg,22), 453 & shp(nshl,maxsh), fresli(npro,22), 454 & fresl(npro,nshl,22), strnrmi(npro), 455 & xnutl(npro), fwr(npro), 456 & xmij(npro,6), xlij(npro,6), 457 & xdenli(npro,ngauss), xnumli(npro,ngauss), 458 & xden(idim), xnum(idim), 459 & lfathe(npro,maxsh) 460 461 call local (fres, fresl, ien, 22, 'gather ') 462 463 do intp = 1, ngauss 464 465 fresli = zero 466 do i = 1, nenl 467 do j = 1, 22 468 fresli(:,j) = fresli(:,j) + shp(i,intp)*fresl(:,i,j) 469 enddo 470 enddo 471 472 strnrmi(:) = sqrt( 473 & two * (fresli(:,10)**2 + fresli(:,11)**2 + fresli(:,12)**2) 474 & + four * ( fresli(:,13)**2 + fresli(:,14)**2 + 475 & fresli(:,15)**2 ) ) 476 477 fwr = fwr1 * fresli(:,22) * strnrmi(:) 478 479 xmij(:,1) = -fwr 480 & * fresli(:,10) + fresli(:,16) 481 482 xmij(:,2) = -fwr 483 & * fresli(:,11) + fresli(:,17) 484 485 xmij(:,3) = -fwr 486 & * fresli(:,12) + fresli(:,18) 487 488 xmij(:,4) = -fwr * fresli(:,13) + fresli(:,19) 489 xmij(:,5) = -fwr * fresli(:,14) + fresli(:,20) 490 xmij(:,6) = -fwr * fresli(:,15) + fresli(:,21) 491 492 fresli(:,22) = one/fresli(:,22) 493 494 xlij(:,1) = fresli(:,4) - fresli(:,1) * 495 & fresli(:,1) * fresli(:,22) 496 xlij(:,2) = fresli(:,5) - fresli(:,2) * 497 & fresli(:,2) * fresli(:,22) 498 xlij(:,3) = fresli(:,6) - fresli(:,3) * 499 & fresli(:,3) * fresli(:,22) 500 xlij(:,4) = fresli(:,7) - fresli(:,1) * 501 & fresli(:,2) * fresli(:,22) 502 xlij(:,5) = fresli(:,8) - fresli(:,1) * 503 & fresli(:,3) * fresli(:,22) 504 xlij(:,6) = fresli(:,9) - fresli(:,2) * 505 & fresli(:,3) * fresli(:,22) 506 507 xnumli(:,intp) = xlij(:,1) * xmij(:,1) 508 & + xlij(:,2) * xmij(:,2) 509 & + xlij(:,3) * xmij(:,3) 510 & + two * (xlij(:,4) * xmij(:,4) 511 & + xlij(:,5) * xmij(:,5) 512 & + xlij(:,6) * xmij(:,6)) 513 xdenli(:,intp) = xmij(:,1) * xmij(:,1) 514 & + xmij(:,2) * xmij(:,2) 515 & + xmij(:,3) * xmij(:,3) 516 & + two * (xmij(:,4) * xmij(:,4) 517 & + xmij(:,5) * xmij(:,5) 518 & + xmij(:,6) * xmij(:,6)) 519 520 521 enddo ! End loop over quad pts 522 523c... For debugging 524 525c do nel = 1, npro 526c do intp = 1, ngauss 527 528c if ( mod(lfathe(nel,intp),2) .eq. 0 )then 529c xnumli(nel,intp) = 2.0 530c xdenli(nel,intp) = 3.0 531c else 532c xnumli(nel,intp) = 5.0 533c xdenli(nel,intp) = 2.5 534c endif 535 536c enddo 537c enddo 538 539c... End of debugging stuff. 540 541 542 do intp = 1, ngauss 543 do nel = 1, npro 544 545 xden( lfathe(nel,intp) ) = xden( lfathe(nel,intp) ) + 546 & xdenli(nel,intp) 547 548 xnum( lfathe(nel,intp) ) = xnum( lfathe(nel,intp) ) + 549 & xnumli(nel,intp) 550 551 enddo 552 enddo 553 554 return 555 556 end 557 558 559 560 561 562 563 564 565 566