1 subroutine getdmc (y, shgl, shp, 2 & iper, ilwork, 3 & nsons, ifath, x) 4 5 use pointer_data 6 7 use quadfilt ! This module gives us shglf(maxtp,nsd,maxsh,ngaussf), 8c shpf(maxtp,maxsh,ngaussf), and Qwtf(maxtp,ngaussf). 9c Shpf and shglf are the shape funciotns and their 10c gradient evaluated using the quadrature rule desired 11c for computing the dmod. Qwt contains the weights of the 12c quad. points. 13 14 15 16 include "common.h" 17 include "mpif.h" 18 include "auxmpi.h" 19 20c 21 dimension fres(nshg,24), fwr(nshg), 22 & strnrm(nshg), cdelsq(nshg), 23 & xnum(nshg), xden(nshg), 24 & xmij(nshg,6), xlij(nshg,6), 25 & xnude(nfath,2), xnuder(nfath,2), 26 & nsons(nfath), 27 & strl(numel,maxnint), 28 & y(nshg,5), 29 & ifath(nshg), iper(nshg), 30 & ilwork(nlwork),! xmudmi(numel,ngauss), 31 & x(numnp,3), 32 & shgl(MAXTOP,nsd,maxsh,MAXQPT), shp(MAXTOP,maxsh,MAXQPT) 33c$$$ & ,xnutf(nfath) must be uncommmented for diags at bottom 34c 35c 36c setup the weights for time averaging of cdelsq (now in quadfilt module) 37c 38 denom=max(1.0d0*(lstep),one) 39 if(dtavei.lt.0) then 40 wcur=one/denom 41 else 42 wcur=dtavei 43 endif 44 whist=1.0-wcur 45c 46c hack in an interesting velocity field (uncomment to test dmod) 47c 48c y(:,5) = 1.0 49c y(:,2) = 2.0*x(:,1) - 3*x(:,2) 50c y(:,3) = 3.0*x(:,1) + 4.0*x(:,2) 51c y(:,4) = 4.0*x(:,1) + x(:,2) + x(:,3) 52c y(:,1) = Rgas * y(:,5) ! Necessary to make model suitable suitable 53 ! for the incompressible case. 54c 55 56 57 fres = zero 58c 59c y(:,5) = 1.0 60c y(:,1) = 2.0*x(:,1) - 3*x(:,2) 61c y(:,2) = 3.0*x(:,1) + 4.0*x(:,2) 62c y(:,3) = 4.0*x(:,1) + x(:,2) + x(:,3) 63c y(:,4) = 1.0 ! Necessary to make model suitable suitable 64 ! for the incompressible case. 65c 66 67 68 intrul=intg(1,itseq) 69 intind=intpt(intrul) 70 71 do iblk = 1,nelblk 72 lcsyst = lcblk(3,iblk) 73 iel = lcblk(1,iblk) !Element number where this block begins 74 npro = lcblk(1,iblk+1) - iel 75 lelCat = lcblk(2,iblk) 76 nenl = lcblk(5,iblk) 77 nshl = lcblk(10,iblk) 78 inum = iel + npro - 1 79 80 ngauss = nint(lcsyst) 81 ngaussf = nintf(lcsyst) 82 83 call asithf (y, x, strl(iel:inum,:), mien(iblk)%p, fres, 84 & shglf(lcsyst,:,1:nshl,:), 85 & shpf(lcsyst,1:nshl,:),Qwtf(lcsyst,1:ngaussf)) 86 87 enddo 88c 89 90 do iblk = 1,nelblk 91 lcsyst = lcblk(3,iblk) 92 iel = lcblk(1,iblk) !Element number where this block begins 93 npro = lcblk(1,iblk+1) - iel 94 lelCat = lcblk(2,iblk) 95 nenl = lcblk(5,iblk) 96 nshl = lcblk(10,iblk) 97 inum = iel + npro - 1 98 99 ngauss = nint(lcsyst) 100 ngaussf = nintf(lcsyst) 101 102 if (ngaussf .ne. ngauss) then 103 call getstrl (y, x, mien(iblk)%p, 104 & strl(iel:inum,:), shgl(lcsyst,:,1:nshl,:), 105 & shp(lcsyst,1:nshl,:)) 106 endif 107 108 enddo 109c 110c 111C must fix for abc and dynamic model 112c if(iabc==1) !are there any axisym bc's 113c & call rotabc(res, iBC, 'in ') 114c 115 if(numpe>1) call commu (fres, ilwork, 24, 'in ') 116!proc-masters = proc-masters + proc-slaves 117c 118c account for periodicity in filtered variables 119c 120 do j = 1,nshg 121 i = iper(j) 122 if (i .ne. j) then 123 fres(i,:) = fres(i,:) + fres(j,:) ! masters = masters + slaves 124 endif 125 enddo 126 do j = 1,nshg 127 i = iper(j) 128 if (i .ne. j) then 129 fres(j,:) = fres(i,:) !slaves get copy of the complete master 130 endif 131 enddo 132 133 if (nfath .eq. 0) then 134 if(numpe>1) call commu (fres, ilwork, 24, 'out') 135c again we will need a call 136c if(iabc==1) !are there any axisym bc's 137c & rotabc(y, iBC, 'out') 138 endif 139 140 fres(:,23) = one / fres(:,23) 141 do j = 1,22 142 fres(:,j) = fres(:,j) * fres(:,23) !"solve" 143 enddo 144c fres(:,24) = fres(:,24) * fres(:,23) 145c 146c.....at this point fres is really all of our filtered quantities 147c at the nodes 148c 149 150 strnrm = sqrt( 151 & two * (fres(:,10)**2 + fres(:,11)**2 + fres(:,12)**2) 152 & + four * ( fres(:,13)**2 + fres(:,14)**2 + fres(:,15)**2 ) ) 153 154 fwr = fwr1 * fres(:,22) * strnrm 155 156 157 158c... There is a difference in xmij(:,1:3) between the compressible 159c... and incompressible forms of the model. 160 161 162 if(matflg(1,1).eq.0) then ! compressible 163 164 xmij(:,1) = -fwr 165 & * pt33 * (two*fres(:,10) - fres(:,11) - fres(:,12)) 166 & + pt33 * (two*fres(:,16) - fres(:,17) - fres(:,18)) 167 xmij(:,2) = -fwr 168 & * pt33 * (two*fres(:,11) - fres(:,10) - fres(:,12)) 169 & + pt33 * (two*fres(:,17) - fres(:,16) - fres(:,18)) 170 xmij(:,3) = -fwr 171 & * pt33 * (two*fres(:,12) - fres(:,10) - fres(:,11)) 172 & + pt33 * (two*fres(:,18) - fres(:,16) - fres(:,17)) 173 174 else 175 176 xmij(:,1) = -fwr 177 & * fres(:,10) + fres(:,16) 178 xmij(:,2) = -fwr 179 & * fres(:,11) + fres(:,17) 180 xmij(:,3) = -fwr 181 & * fres(:,12) + fres(:,18) 182 183 endif 184 185 xmij(:,4) = -fwr * fres(:,13) + fres(:,19) 186 xmij(:,5) = -fwr * fres(:,14) + fres(:,20) 187 xmij(:,6) = -fwr * fres(:,15) + fres(:,21) 188 189 fres(:,22) = one / fres(:,22) 190 191 xlij(:,1) = fres(:,4) - fres(:,1) * fres(:,1) * fres(:,22) 192 xlij(:,2) = fres(:,5) - fres(:,2) * fres(:,2) * fres(:,22) 193 xlij(:,3) = fres(:,6) - fres(:,3) * fres(:,3) * fres(:,22) 194 xlij(:,4) = fres(:,7) - fres(:,1) * fres(:,2) * fres(:,22) 195 xlij(:,5) = fres(:,8) - fres(:,1) * fres(:,3) * fres(:,22) 196 xlij(:,6) = fres(:,9) - fres(:,2) * fres(:,3) * fres(:,22) 197 198 xnum = xlij(:,1) * xmij(:,1) + xlij(:,2) * xmij(:,2) 199 & + xlij(:,3) * xmij(:,3) 200 & + two * (xlij(:,4) * xmij(:,4) + xlij(:,5) * xmij(:,5) 201 & + xlij(:,6) * xmij(:,6)) 202 xden = xmij(:,1) * xmij(:,1) + xmij(:,2) * xmij(:,2) 203 & + xmij(:,3) * xmij(:,3) 204 & + two * (xmij(:,4) * xmij(:,4) + xmij(:,5) * xmij(:,5) 205 & + xmij(:,6) * xmij(:,6)) 206 xden = two * xden 207 208c zero on processor periodic nodes so that they will not be added twice 209 do j = 1,numnp 210 i = iper(j) 211 if (i .ne. j) then 212 xnum(j) = zero 213 xden(j) = zero 214 endif 215 enddo 216 217 if (numpe.gt.1 .and. nsons(1).gt.1) then 218 219 numtask = ilwork(1) 220 itkbeg = 1 221 222c zero the nodes that are "solved" on the other processors 223 do itask = 1, numtask 224 225 iacc = ilwork (itkbeg + 2) 226 numseg = ilwork (itkbeg + 4) 227 228 if (iacc .eq. 0) then 229 do is = 1,numseg 230 isgbeg = ilwork (itkbeg + 3 + 2*is) 231 lenseg = ilwork (itkbeg + 4 + 2*is) 232 isgend = isgbeg + lenseg - 1 233 xnum(isgbeg:isgend) = zero 234 xden(isgbeg:isgend) = zero 235 enddo 236 endif 237 238 itkbeg = itkbeg + 4 + 2*numseg 239 240 enddo 241 242 endif 243c 244c Description of arrays. Each processor has an array of length equal 245c to the total number of fathers times 2 xnude(nfathers,2). One to collect 246c the numerator and one to collect the denominator. There is also an array 247c of length nshg on each processor which tells the father number of each 248c on processor node, ifath(nnshg). Finally, there is an arry of length 249c nfathers to tell the total (on all processors combined) number of sons 250c for each father. 251c 252c Now loop over nodes and accumlate the numerator and the denominator 253c to the father nodes. Only on processor addition at this point. 254c Note that serrogate fathers are collect some for the case where some 255c sons are on another processor 256c 257 ihomog=1 258 if(maxval(nsons).eq.1) ihomog=0 259 if(ihomog.eq.1) then 260 xnude = zero 261 do i = 1,nshg 262 xnude(ifath(i),1) = xnude(ifath(i),1) + xnum(i) 263 xnude(ifath(i),2) = xnude(ifath(i),2) + xden(i) 264 enddo 265 else ! no homogeneity assumed 266 xnude(:,1)=xnum(:) 267 xnude(:,2)=xden(:) 268 endif 269 270c 271c Now the true fathers and serrogates combine results and update 272c each other. 273c 274c ONLY DO THIS IF 1) multiprocessors AND 2) homogenous directions exist 275 276 if(numpe .gt. 1 .and. ihomog.eq.1 )then 277 call drvAllreduce(xnude, xnuder,2*nfath) 278c 279c xnude is the sum of the sons for each father on this processor 280c 281c xnuder is the sum of the sons for each father on all processor combined 282c (the same as if we had not partitioned the mesh for each processor) 283c 284c For each father we have precomputed the number of sons (including 285c the sons off processor). 286c 287c Now divide by number of sons to get the average (not really necessary 288c for dynamic model since ratio will cancel nsons at each father) 289c 290c xnuder(:,1) = xnuder(:,1) ! / nsons(:) 291c xnuder(:,2) = xnuder(:,2) ! / nsons(:) 292c 293c the next line is c \Delta^2 294c 295 numNden(:,1) = whist*numNden(:,1)+wcur*xnuder(ifath(:),1) 296 numNden(:,2) = whist*numNden(:,2)+wcur*xnuder(ifath(:),2) 297 cdelsq(:) = numNden(:,1) / (numNden(:,2) + 1.d-09) 298c note that we have whist and wcur in here to allow for both time 299c averaging to be used in conjunction with spatial homogenous averaging 300 301 else 302c 303c the next line is c \Delta^2, not nu_T but we want to save the 304c memory 305c 306 307c$$$ write(400+myrank,555) (numNden(j*500,1),j=1,5) 308c$$$ write(410+myrank,555) (numNden(j*500,2),j=1,5) 309c$$$ write(500+myrank,555) (numNden(j+500,1),j=1,5) 310c$$$ write(510+myrank,555) (numNden(j+500,2),j=1,5) 311c$$$ 312c$$$ write(430+myrank,555) (xnude(j*500,1),j=1,5) 313c$$$ write(440+myrank,555) (xnude(j*500,2),j=1,5) 314c$$$ write(530+myrank,555) (xnude(j+500,1),j=1,5) 315c$$$ write(540+myrank,555) (xnude(j+500,2),j=1,5) 316 317 numNden(:,1) = whist*numNden(:,1)+wcur*xnude(ifath(:),1) 318 numNden(:,2) = whist*numNden(:,2)+wcur*xnude(ifath(:),2) 319 cdelsq(:) = numNden(:,1) / (numNden(:,2) + 1.d-09) 320 321c 322c to get started we hold cdelsq fixed 323c 324c cdelsq(:) = 2.27e-4 325c cdelsq(:) = 0 326 327 328c$$$ 329c$$$ write(450+myrank,555) (cdelsq(j*500),j=1,5) 330c$$$ write(460+myrank,555) (y(j*500,1),j=1,5) 331c$$$ write(470+myrank,555) (y(j*500,2),j=1,5) 332c$$$ write(480+myrank,555) (y(j*500,3),j=1,5) 333c$$$ write(490+myrank,555) (strnrm(j*500),j=1,5) 334c$$$ 335c$$$ write(550+myrank,555) (cdelsq(j+500),j=1,5) 336c$$$ write(560+myrank,555) (y(j+500,1),j=1,5) 337c$$$ write(570+myrank,555) (y(j+500,2),j=1,5) 338c$$$ write(580+myrank,555) (y(j+500,3),j=1,5) 339c$$$ write(590+myrank,555) (strnrm(j+500),j=1,5) 340c$$$ 341c$$$ call flush(400+myrank) 342c$$$ call flush(410+myrank) 343c$$$ call flush(430+myrank) 344c$$$ call flush(440+myrank) 345c$$$ call flush(450+myrank) 346c$$$ call flush(460+myrank) 347c$$$ call flush(470+myrank) 348c$$$ call flush(480+myrank) 349c$$$ call flush(490+myrank) 350c$$$ call flush(500+myrank) 351c$$$ call flush(510+myrank) 352c$$$ call flush(530+myrank) 353c$$$ call flush(540+myrank) 354c$$$ call flush(550+myrank) 355c$$$ call flush(560+myrank) 356c$$$ call flush(570+myrank) 357c$$$ call flush(580+myrank) 358c$$$ call flush(590+myrank) 359 endif 360 555 format(5(2x,e14.7)) 361 362c $$$$$$$$$$$$$$$$$$$$$$$$$$$ 363 tmp1 = MINVAL(cdelsq) 364 tmp2 = MAXVAL(cdelsq) 365 if(numpe>1) then 366 call MPI_REDUCE (tmp1, tmp3, 1,MPI_DOUBLE_PRECISION, 367 & MPI_MIN, master, MPI_COMM_WORLD, ierr) 368 call MPI_REDUCE (tmp2, tmp4, 1, MPI_DOUBLE_PRECISION, 369 & MPI_MAX, master, MPI_COMM_WORLD, ierr) 370 tmp1=tmp3 371 tmp2=tmp4 372 endif 373 if (myrank .EQ. master) then !print CDelta^2 range 374 write(34,*)lstep,tmp1,tmp2 375 call flush(34) 376 endif 377c $$$$$$$$$$$$$$$$$$$$$$$$$$$ 378 379 if (myrank .eq. master) then 380 write(*,*)'xnut=',sum(cdelsq)/nshg 381 write(*,*) 'cdelsq=', cdelsq(1),cdelsq(2) 382 endif 383 384 do iblk = 1,nelblk 385 lcsyst = lcblk(3,iblk) 386 iel = lcblk(1,iblk) 387 npro = lcblk(1,iblk+1) - iel 388 lelCat = lcblk(2,iblk) 389 inum = iel + npro - 1 390 391 ngauss = nint(lcsyst) 392 393 call scatnu (mien(iblk)%p, strl(iel:inum,:), 394 & mxmudmi(iblk)%p,cdelsq,shp(lcsyst,1:nshl,:)) 395 enddo 396c $$$$$$$$$$$$$$$$$$$$$$$$$$$ 397c$$$ tmp1 = MINVAL(xmudmi) 398c$$$ tmp2 = MAXVAL(xmudmi) 399c$$$ if(numpe>1) then 400c$$$ call MPI_REDUCE (tmp1, tmp3, 1, MPI_DOUBLE_PRECISION, 401c$$$ & MPI_MIN, master, MPI_COMM_WORLD, ierr) 402c$$$ call MPI_REDUCE (tmp2, tmp4, 1, MPI_DOUBLE_PRECISION, 403c$$$ & MPI_MAX, master, MPI_COMM_WORLD, ierr) 404c$$$ tmp1=tmp3 405c$$$ tmp2=tmp4 406c$$$ endif 407c$$$ if (myrank .EQ. master) then 408c$$$ write(35,*) lstep,tmp1,tmp2 409c$$$ call flush(35) 410c$$$ endif 411c $$$$$$$$$$$$$$$$$$$$$$$$$$$ 412 413c$$$c 414c$$$c if flag set, write a restart file with info (reuse xmij's memory) 415c$$$c 416c$$$ if(irs.eq.11) then 417c$$$ lstep=999 418c$$$ xmij(:,1)=xnum(:) 419c$$$ xmij(:,2)=xden(:) 420c$$$ xmij(:,3)=cdelsq(:) 421c$$$ xmij(:,5)=xlij(:,4) !leave M_{12} in 4 and put L_{12} here 422c$$$ call restar('out ',xmij,xlij) !also dump all of L_{ij} in ac 423c$$$ stop 424c$$$ endif 425c$$$ if(lhs.eq.1) then 426c$$$ lstepsafe=lstep 427c$$$ lstep=999 428c$$$ 429c$$$ xmij(:,1)=xnum(:) 430c$$$ xmij(:,2)=xden(:) 431c$$$ xmij(:,3)=cdelsq(:) 432c$$$ xmij(:,4)=xmij(:,6) !put M_{23} in 4 433c$$$ xmij(:,5)=xlij(:,6) !put L_{23} here 434c$$$ call restar('out ',xmij,xlij) !also dump all of L_{ij} in ac 435c$$$ lstep=lstepsafe 436c$$$ endif 437c$$$c 438c$$$c local clipping moved to scatnu with the creation of mxmudmi pointers 439c$$$c 440c$$$c$$$ rmu=datmat(1,2,1) 441c$$$c$$$ xmudmi=min(xmudmi,1000.0*rmu) !don't let it get larger than 1000 mu 442c$$$c$$$ xmudmi=max(xmudmi, -rmu) ! don't let (xmudmi + mu) < 0 443c$$$c stop !uncomment to test dmod 444c$$$c 445c$$$ 446c$$$ 447c$$$c write out the nodal values of xnut (estimate since we don't calc strain 448c$$$c there and must use the filtered strain). 449c$$$c 450c$$$ 451c$$$ if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then 452c$$$c 453c$$$c collect the average strain into xnude(2) 454c$$$c 455c$$$ xnude(:,2) = zero 456c$$$ do i = 1,numnp 457c$$$ xnude(ifath(i),2) = xnude(ifath(i),2) + strnrm(i) 458c$$$ enddo 459c$$$ 460c$$$ if(numpe .gt. 1 .and. ihomog.eq.1 )then 461c$$$ call drvAllreduce(xnude(:,2), xnuder(:,2),nfath) 462c$$$ else 463c$$$ xnuder=xnude 464c$$$ endif 465c$$$c 466c$$$c nut= cdelsq * |S| 467c$$$c 468c$$$ xnutf=xnuder(:,1)*xnuder(:,2)/nsons(:) ! off by one 469c$$$c 470c$$$c collect the x and y coords into xnude 471c$$$c 472c$$$ xnude = zero 473c$$$ do i = 1,numnp 474c$$$ xnude(ifath(i),1) = xnude(ifath(i),1) + x(i,3) 475c$$$ xnude(ifath(i),2) = xnude(ifath(i),2) + x(i,2) 476c$$$ enddo 477c$$$ 478c$$$ if(numpe .gt. 1 .and. nsons(1).ne.1 )then 479c$$$ call drvAllreduce(xnude, xnuder,2*nfath) 480c$$$ xnuder(:,1)=xnuder(:,1)/nsons(:) 481c$$$ xnuder(:,2)=xnuder(:,2)/nsons(:) 482c$$$ else 483c$$$ xnuder=xnude 484c$$$ endif 485c$$$c 486c$$$c xnude is the sum of the sons for each father on this processor 487c$$$c 488c$$$ if((myrank.eq.master)) then 489c$$$ do i=1,nfath ! cdelsq * |S| 490c$$$ write(444,*) xnuder(i,1),xnuder(i,2),xnutf(i) 491c$$$ enddo 492c$$$ call flush(444) 493c$$$ endif 494c$$$ endif 495 496 return 497 end 498