1 subroutine itrdrv (y, ac, 2 & uold, x, 3 & iBC, BC, 4 & iper, ilwork, shp, 5 & shgl, shpb, shglb, 6 & ifath, velbar, nsons ) 7c 8c---------------------------------------------------------------------- 9c 10c This iterative driver is the semi-discrete, predictor multi-corrector 11c algorithm. It contains the Hulbert Generalized Alpha method which 12c is 2nd order accurate for Rho_inf from 0 to 1. The method can be 13c made first-order accurate by setting Rho_inf=-1. It uses CGP and 14c GMRES iterative solvers. 15c 16c working arrays: 17c y (nshg,ndof) : Y variables 18c x (nshg,nsd) : node coordinates 19c iBC (nshg) : BC codes 20c BC (nshg,ndofBC) : BC constraint parameters 21c iper (nshg) : periodicity table 22c 23c 24c Zdenek Johan, Winter 1991. (Fortran 90) 25c Alberto Figueroa, Winter 2004. CMM-FSI 26c Irene Vignon, Fall 2004. Impedance BC 27c---------------------------------------------------------------------- 28c 29 use pvsQbi !gives us splag (the spmass at the end of this run 30 use specialBC !gives us itvn 31 use timedata !allows collection of time series 32 use convolImpFlow !for Imp bc 33 use convolRCRFlow !for RCR bc 34!MR CHANGE 35 use turbsa ! used to access d2wall 36!MR CHANGE END 37 38c use readarrays !reads in uold and acold 39 40 include "common.h" 41 include "mpif.h" 42 include "auxmpi.h" 43c 44 45 46 real*8 y(nshg,ndof), ac(nshg,ndof), 47 & yold(nshg,ndof), acold(nshg,ndof), 48 & u(nshg,nsd), uold(nshg,nsd), 49 & x(numnp,nsd), solinc(nshg,ndof), 50 & BC(nshg,ndofBC), tf(nshg,ndof), 51 & GradV(nshg,nsdsq) 52 53c 54 real*8 res(nshg,ndof) 55c 56 real*8 shp(MAXTOP,maxsh,MAXQPT), 57 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 58 & shpb(MAXTOP,maxsh,MAXQPT), 59 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 60c 61 integer rowp(nshg,nnz), colm(nshg+1), 62 & iBC(nshg), 63 & ilwork(nlwork), 64 & iper(nshg), ifuncs(6) 65 66 real*8 vbc_prof(nshg,3) 67 68 integer stopjob 69 character*10 cname2 70 character*5 cname 71c 72c stuff for dynamic model s.w.avg and wall model 73c 74 dimension ifath(numnp), velbar(nfath,ndof), nsons(nfath) 75 76 dimension wallubar(2),walltot(2) 77c 78c.... For Farzin's Library 79c 80 integer eqnType, prjFlag, presPrjFlag, verbose 81c 82 real*8, allocatable, dimension(:,:) :: aperm, atemp, atempS 83 real*8, allocatable, dimension(:,:,:) :: apermS 84 85 real*8, allocatable, dimension(:,:) :: lhsP, lhsK, lhsS 86 real*8 almit, alfit, gamit 87c 88 character*1024 servername 89 character*20 fname1,fmt1 90 character*20 fname2,fmt2 91 character*60 fnamepold, fvarts 92 character*4 fname4c ! 4 characters 93 integer iarray(50) ! integers for headers 94 integer isgn(ndof), isgng(ndof) 95 96!MR CHANGE 97! real*8 rerr(nshg,10), ybar(nshg,13) ! including 7 sq. terms with 3 cross terms of uv, uw and vw 98! real*8 rerr(nshg,10), ybar(nshg,12) ! including 7 sq. terms with 3 cross terms of uv, uw and vw 99 real*8 rerr(nshg,10) 100 real*8, allocatable, dimension(:,:) :: ybar, strain, vorticity 101 real*8, allocatable, dimension(:,:) :: wallssVec, wallssVecbar 102!MR CHANGE 103 104 real*8 tcorecp(2), tcorecpscal(2) 105 106 integer, allocatable, dimension(:) :: ivarts 107 integer, allocatable, dimension(:) :: ivartsg 108 integer, allocatable, dimension(:) :: iv_rank 109 real*8, allocatable, dimension(:) :: vartssoln 110 real*8, allocatable, dimension(:) :: vartssolng 111 real*8, allocatable, dimension(:,:,:) :: yphbar 112 real*8 CFLworst(numel) 113 114!MR CHANGE 115 integer :: iv_rankpernode, iv_totnodes, iv_totcores 116 integer :: iv_node, iv_core, iv_thread 117!MR CHANGE 118 119 impistat = 0 120 impistat2 = 0 121 iISend = 0 122 iISendScal = 0 123 iIRecv = 0 124 iIRecvScal = 0 125 iWaitAll = 0 126 iWaitAllScal = 0 127 iAllR = 0 128 iAllRScal = 0 129 rISend = zero 130 rISendScal = zero 131 rIRecv = zero 132 rIRecvScal = zero 133 rWaitAll = zero 134 rWaitAllScal = zero 135 rAllR = zero 136 rAllRScal = zero 137 rCommu = zero 138 rCommuScal = zero 139 140 call initmemstat() 141 142c 143c hack SA variable 144c 145cHack BC(:,7)=BC(:,7)*0.001 146cHack if(lstep.eq.0) y(:,6)=y(:,6)*0.001 147 call SolverLicenseServer(servername) 148c 149c only master should be verbose 150c 151 152 if(numpe.gt.0 .and. myrank.ne.master)iverbose=0 153c 154 155 lskeep=lstep 156 157 inquire(file='xyzts.dat',exist=exts) 158 if(exts) then 159 160 open(unit=626,file='xyzts.dat',status='old') 161 read(626,*) ntspts, freq, tolpt, iterat, varcod 162 call sTD ! sets data structures 163 164 do jj=1,ntspts ! read coordinate data where solution desired 165 read(626,*) ptts(jj,1),ptts(jj,2),ptts(jj,3) 166 enddo 167 close(626) 168 169 statptts(:,:) = 0 170 parptts(:,:) = zero 171 varts(:,:) = zero 172 173 allocate (ivarts(ntspts*ndof)) 174 allocate (ivartsg(ntspts*ndof)) 175 allocate (iv_rank(ntspts)) 176 allocate (vartssoln(ntspts*ndof)) 177 allocate (vartssolng(ntspts*ndof)) 178 179 iv_rankpernode = iv_rankpercore*iv_corepernode 180 iv_totnodes = numpe/iv_rankpernode 181 iv_totcores = iv_corepernode*iv_totnodes 182 if (myrank .eq. 0) then 183 write(*,*) 'Info for probes:' 184 write(*,*) ' Ranks per core:',iv_rankpercore 185 write(*,*) ' Cores per node:',iv_corepernode 186 write(*,*) ' Ranks per node:',iv_rankpernode 187 write(*,*) ' Total number of nodes:',iv_totnodes 188 write(*,*) ' Total number of cores',iv_totcores 189 endif 190 191! if (myrank .eq. numpe-1) then 192 do jj=1,ntspts 193 194 ! Compute the adequate rank which will take care of probe jj 195 jjm1 = jj-1 196 iv_node = (iv_totnodes-1)-mod(jjm1,iv_totnodes) 197 iv_core = (iv_corepernode-1) - mod((jjm1 - 198 & mod(jjm1,iv_totnodes))/iv_totnodes,iv_corepernode) 199 iv_thread = (iv_rankpercore-1) - mod((jjm1- 200 & (mod(jjm1,iv_totcores)))/iv_totcores,iv_rankpercore) 201 iv_rank(jj) = iv_node*iv_rankpernode 202 & + iv_core*iv_rankpercore 203 & + iv_thread 204 205 if(myrank == 0) then 206 write(*,*) ' Probe', jj, 'handled by rank', 207 & iv_rank(jj), ' on node', iv_node 208 endif 209 210 ! Verification just in case 211 if(iv_rank(jj) .lt.0 .or. iv_rank(jj) .ge. numpe) then 212 write(*,*) 'WARNING: iv_rank(',jj,') is ', iv_rank(jj), 213 & ' and reset to numpe-1' 214 iv_rank(jj) = numpe-1 215 endif 216 217 ! Open the varts files 218 if(myrank == iv_rank(jj)) then 219 fvarts='varts/varts' 220 fvarts=trim(fvarts)//trim(cname2(jj)) 221 fvarts=trim(fvarts)//trim(cname2(lstep)) 222 fvarts=trim(fvarts)//'.dat' 223 fvarts=trim(fvarts) 224 open(unit=1000+jj, file=fvarts, status='unknown') 225 endif 226 enddo 227! endif 228 229 endif 230c 231c.... open history and aerodynamic forces files 232c 233 if (myrank .eq. master) then 234! open (unit=ihist, file=fhist, status='unknown') 235 open (unit=iforce, file=fforce, status='unknown') 236 open (unit=76, file="fort.76", status='unknown') 237 if(numImpSrfs.gt.0 .or. numRCRSrfs.gt.0) then 238 fnamepold = 'pold' 239 fnamepold = trim(fnamepold)//trim(cname2(lstep)) 240 fnamepold = trim(fnamepold)//'.dat' 241 fnamepold = trim(fnamepold) 242 open (unit=8177, file=fnamepold, status='unknown') 243 endif 244 endif 245c 246c.... initialize 247c 248 ifuncs(:) = 0 ! func. evaluation counter 249 istep = 0 250 yold = y 251 acold = ac 252 253 rerr = zero 254 255 if(ierrcalc.eq.1 .or. ioybar.eq.1) then ! we need ybar for error too 256 if (ivort == 1) then 257 allocate(ybar(nshg,17)) ! more space for vorticity if requested 258 else 259 allocate(ybar(nshg,13)) 260 endif 261 ybar = zero ! Initialize ybar to zero, which is essential 262 endif 263 264 if(ivort == 1) then 265 allocate(strain(nshg,6)) 266 allocate(vorticity(nshg,5)) 267 endif 268 269 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 270 allocate(wallssVec(nshg,3)) 271 if (ioybar .eq. 1) then 272 allocate(wallssVecbar(nshg,3)) 273 wallssVecbar = zero ! Initialization important if mean wss computed 274 endif 275 endif 276 277! both nstepsincycle and nphasesincycle needs to be set 278 if(nstepsincycle.eq.0) nphasesincycle = 0 279 if(nphasesincycle.ne.0) then 280! & allocate(yphbar(nshg,5,nphasesincycle)) 281 if (ivort == 1) then 282 allocate(yphbar(nshg,15,nphasesincycle)) ! more space for vorticity 283 else 284 allocate(yphbar(nshg,11,nphasesincycle)) 285 endif 286 yphbar = zero 287 endif 288 289!MR CHANGE END 290 291 vbc_prof(:,1:3) = BC(:,3:5) 292 if(iramp.eq.1) then 293 call BCprofileInit(vbc_prof,x) 294 endif 295 296c 297c.... ---------------> initialize Farzin's Library <--------------- 298c 299c.... assign parameter values 300c 301 do i = 1, 100 302 numeqns(i) = i 303 enddo 304 nKvecs = Kspace 305 prjFlag = iprjFlag 306 presPrjFlag = ipresPrjFlag 307 verbose = iverbose 308c 309c.... determine how many scalar equations we are going to need to solve 310c 311 nsolt=mod(impl(1),2) ! 1 if solving temperature 312 nsclrsol=nsolt+nsclr ! total number of scalars solved At 313 ! some point we probably want to create 314 ! a map, considering stepseq(), to find 315 ! what is actually solved and only 316 ! dimension lhs to the appropriate 317 ! size. (see 1.6.1 and earlier for a 318 ! "failed" attempt at this). 319 320 321 nsolflow=mod(impl(1),100)/10 ! 1 if solving flow 322 323c 324c.... Now, call Farzin's lesNew routine to initialize 325c memory space 326c 327 call genadj(colm, rowp, icnt ) ! preprocess the adjacency list 328 329 nnz_tot=icnt ! this is exactly the number of non-zero blocks on 330 ! this proc 331 332 if (nsolflow.eq.1) then 333 lesId = numeqns(1) 334 eqnType = 1 335 nDofs = 4 336 call myfLesNew( lesId, 41994, 337 & eqnType, 338 & nDofs, minIters, maxIters, 339 & nKvecs, prjFlag, nPrjs, 340 & presPrjFlag, nPresPrjs, epstol(1), 341 & prestol, verbose, statsflow, 342 & nPermDims, nTmpDims, servername ) 343 344 allocate (aperm(nshg,nPermDims)) 345 allocate (atemp(nshg,nTmpDims)) 346 allocate (lhsP(4,nnz_tot)) 347 allocate (lhsK(9,nnz_tot)) 348 349 call readLesRestart( lesId, aperm, nshg, myrank, lstep, 350 & nPermDims ) 351 352 else 353 nPermDims = 0 354 nTempDims = 0 355 endif 356 357 358 if(nsclrsol.gt.0) then 359 do isolsc=1,nsclrsol 360 lesId = numeqns(isolsc+1) 361 eqnType = 2 362 nDofs = 1 363 presPrjflag = 0 364 nPresPrjs = 0 365 prjFlag = 1 366 indx=isolsc+2-nsolt ! complicated to keep epstol(2) for 367 ! temperature followed by scalars 368 call myfLesNew( lesId, 41994, 369 & eqnType, 370 & nDofs, minIters, maxIters, 371 & nKvecs, prjFlag, nPrjs, 372 & presPrjFlag, nPresPrjs, epstol(indx), 373 & prestol, verbose, statssclr, 374 & nPermDimsS, nTmpDimsS, servername ) 375 enddo 376c 377c Assume all scalars have the same size needs 378c 379 allocate (apermS(nshg,nPermDimsS,nsclrsol)) 380 allocate (atempS(nshg,nTmpDimsS)) !they can all share this 381 allocate (lhsS(nnz_tot,nsclrsol)) 382c 383c actually they could even share with atemp but leave that for later 384c 385 else 386 nPermDimsS = 0 387 nTmpDimsS = 0 388 endif 389c 390c... prepare lumped mass if needed 391c 392 if((flmpr.ne.0).or.(flmpl.ne.0)) call genlmass(x, shp,shgl) 393c 394c.... -----------------> End of initialization <----------------- 395c 396c.....open the necessary files to gather time series 397c 398 lstep0 = lstep+1 399 nsteprcr = nstep(1)+lstep 400c 401c.... loop through the time sequences 402c 403 404 405 do 3000 itsq = 1, ntseq 406 itseq = itsq 407 408CAD tcorecp1 = second(0) 409CAD tcorewc1 = second(-1) 410c 411c.... set up the time integration parameters 412c 413 nstp = nstep(itseq) 414 nitr = niter(itseq) 415 LCtime = loctim(itseq) 416 dtol(:)= deltol(itseq,:) 417 418 call itrSetup ( y, acold ) 419c 420c...initialize the coefficients for the impedance convolution, 421c which are functions of alphaf so need to do it after itrSetup 422 if(numImpSrfs.gt.zero) then 423 call calcImpConvCoef (numImpSrfs, ntimeptpT) 424 endif 425c 426c...initialize the initial condition P(0)-RQ(0)-Pd(0) for RCR BC 427c need ndsurf so should be after initNABI 428 if(numRCRSrfs.gt.zero) then 429 call calcRCRic(y,nsrflistRCR,numRCRSrfs) 430 endif 431c 432c find the last solve of the flow in the step sequence so that we will 433c know when we are at/near end of step 434c 435c ilast=0 436 nitr=0 ! count number of flow solves in a step (# of iterations) 437 do i=1,seqsize 438 if(stepseq(i).eq.0) nitr=nitr+1 439 enddo 440 441 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 442 tcorecp(:) = zero ! used in solfar.f (solflow) 443 tcorecpscal(:) = zero ! used in solfar.f (solflow) 444 if(myrank.eq.0) then 445 tcorecp1 = TMRC() 446 endif 447 448c 449c.... loop through the time steps 450c 451 istop=0 452 rmub=datmat(1,2,1) 453 if(rmutarget.gt.0) then 454 rmue=rmutarget 455 else 456 rmue=datmat(1,2,1) ! keep constant 457 endif 458 459 if(iramp.eq.1) then 460 call BCprofileScale(vbc_prof,BC,yold) ! fix the yold values to the reset BC 461 isclr=1 ! fix scalar 462 do isclr=1,nsclr 463 call itrBCSclr(yold,ac,iBC,BC,iper,ilwork) 464 enddo 465 endif 466 467!MR CHANGE 468! do irank=0,numpe-1 469! if(irank==myrank) then 470! write(*,*) 'NUMTASK - rank,numtask: ',myrank,ilwork(1) 471! endif 472! call MPI_Barrier(MPI_COMM_WORLD,ierr) 473! enddo 474!MR CHANGE END 475 476 do 2000 istp = 1, nstp 477 if(iramp.eq.1) 478 & call BCprofileScale(vbc_prof,BC,yold) 479 480 call rerun_check(stopjob) 481 if(stopjob.ne.0) goto 2001 482 483 xi=istp*1.0/nstp 484 datmat(1,2,1)=rmub*(1.0-xi)+xi*rmue 485c write(*,*) "current mol. visc = ", datmat(1,2,1) 486c.... if we have time varying boundary conditions update the values of BC. 487c these will be for time step n+1 so use lstep+1 488c 489 if(itvn.gt.0) call BCint((lstep+1)*Delt(1), shp, shgl, 490 & shpb, shglb, x, BC, iBC) 491 492c 493c ... calculate the pressure contribution that depends on the history for the Imp. BC 494c 495 if(numImpSrfs.gt.0) then 496 call pHist(poldImp,QHistImp,ImpConvCoef, 497 & ntimeptpT,numImpSrfs) 498 if (myrank.eq.master) 499 & write(8177,*) (poldImp(n), n=1,numImpSrfs) 500 endif 501c 502c ... calc the pressure contribution that depends on the history for the RCR BC 503c 504 if(numRCRSrfs.gt.0) then 505 call CalcHopRCR (Delt(itseq), lstep, numRCRSrfs) 506 call CalcRCRConvCoef(lstep,numRCRSrfs) 507 call pHist(poldRCR,QHistRCR,RCRConvCoef,nsteprcr, 508 & numRCRSrfs) 509 if (myrank.eq.master) 510 & write(8177,*) (poldRCR(n), n=1,numRCRSrfs) 511 endif 512c 513c Decay of scalars 514c 515 if(nsclr.gt.0 .and. tdecay.ne.1) then 516 yold(:,6:ndof)=y(:,6:ndof)*tdecay 517 BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*tdecay 518 endif 519 520 if(nosource.eq.1) BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*0.8 521 522 523 if(iLES.gt.0) then !complicated stuff has moved to 524 !routine below 525 call lesmodels(yold, acold, shgl, shp, 526 & iper, ilwork, rowp, colm, 527 & nsons, ifath, x, 528 & iBC, BC) 529 530 531 endif 532 533c.... set traction BCs for modeled walls 534c 535 if (itwmod.ne.0) then 536 call asbwmod(yold, acold, x, BC, iBC, 537 & iper, ilwork, ifath, velbar) 538 endif 539 540!MR CHANGE 541c 542c.... Determine whether the vorticity field needs to be computed for this time step or not 543c 544 icomputevort = 0 545 if (ivort == 1) then ! Print vorticity = True in solver.inp 546 ! We then compute the vorticity only if we 547 ! 1) we write an intermediate checkpoint 548 ! 2) we reach the last time step and write the last checkpoint 549 ! 3) we accumulate statistics in ybar for every time step 550 ! BEWARE: we need here lstep+1 and istep+1 because the lstep and 551 ! istep gets incremened after the flowsolve, further below 552 if (((irs .ge. 1) .and. (mod(lstep+1, ntout) .eq. 0)) .or. 553 & istep+1.eq.nstep(itseq) .or. ioybar == 1) then 554 icomputevort = 1 555 endif 556 endif 557 558! write(*,*) 'icomputevort: ',icomputevort, ' - istep: ', 559! & istep,' - nstep(itseq):',nstep(itseq),'- lstep:', 560! & lstep, '- ntout:', ntout 561!MR CHANGE 562 563c 564c.... -----------------------> predictor phase <----------------------- 565c 566 call itrPredict(yold, y, acold, ac , uold, u, iBC) 567 call itrBC (y, ac, iBC, BC, iper,ilwork) 568 569 if(nsolt.eq.1) then 570 isclr=0 571 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 572 endif 573 do isclr=1,nsclr 574 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 575 enddo 576 iter=0 577 ilss=0 ! this is a switch thrown on first solve of LS redistance 578 do istepc=1,seqsize 579 icode=stepseq(istepc) 580 if(mod(icode,5).eq.0) then ! this is a solve 581 isolve=icode/10 582 if(icode.eq.0) then ! flow solve (encoded as 0) 583c 584 iter = iter+1 585 ifuncs(1) = ifuncs(1) + 1 586c 587 Force(1) = zero 588 Force(2) = zero 589 Force(3) = zero 590 HFlux = zero 591 lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 592 593 call SolFlow(y, ac, u, 594 & yold, acold, uold, 595 & x, iBC, 596 & BC, res, 597 & nPermDims, nTmpDims, aperm, 598 & atemp, iper, 599 & ilwork, shp, shgl, 600 & shpb, shglb, rowp, 601 & colm, lhsK, lhsP, 602 & solinc, rerr, tcorecp, 603 & GradV) 604 605 else ! scalar type solve 606 if (icode.eq.5) then ! Solve for Temperature 607 ! (encoded as (nsclr+1)*10) 608 isclr=0 609 ifuncs(2) = ifuncs(2) + 1 610 j=1 611 else ! solve a scalar (encoded at isclr*10) 612 isclr=isolve 613 ifuncs(isclr+2) = ifuncs(isclr+2) + 1 614 j=isclr+nsolt 615 if((iLSet.eq.2).and.(ilss.eq.0) 616 & .and.(isclr.eq.2)) then 617 ilss=1 ! throw switch (once per step) 618 y(:,7)=y(:,6) ! redistance field initialized 619 ac(:,7) = zero 620 call itrBCSclr ( y, ac, iBC, BC, iper, 621 & ilwork) 622c 623c....store the flow alpha, gamma parameter values and assigm them the 624c....Backward Euler parameters to solve the second levelset scalar 625c 626 alfit=alfi 627 gamit=gami 628 almit=almi 629 Deltt=Delt(1) 630 Dtglt=Dtgl 631 alfi = 1 632 gami = 1 633 almi = 1 634c Delt(1)= Deltt ! Give a pseudo time step 635 Dtgl = one / Delt(1) 636 endif ! level set eq. 2 637 endif ! deciding between temperature and scalar 638 639 lhs = 1 - min(1,mod(ifuncs(isclr+2)-1, 640 & LHSupd(isclr+2))) 641 642 call SolSclr(y, ac, u, 643 & yold, acold, uold, 644 & x, iBC, 645 & BC, nPermDimsS,nTmpDimsS, 646 & apermS(1,1,j), atempS, iper, 647 & ilwork, shp, shgl, 648 & shpb, shglb, rowp, 649 & colm, lhsS(1,j), 650 & solinc(1,isclr+5), tcorecpscal) 651 652 653 endif ! end of scalar type solve 654 655 else ! this is an update (mod did not equal zero) 656 iupdate=icode/10 ! what to update 657 if(icode.eq.1) then !update flow 658 call itrCorrect ( y, ac, u, solinc, iBC) 659 call itrBC (y, ac, iBC, BC, iper, ilwork) 660 else ! update scalar 661 isclr=iupdate !unless 662 if(icode.eq.6) isclr=0 663 if(iRANS.lt.-100)then ! RANS 664 call itrCorrectSclrPos(y,ac,solinc(1,isclr+5)) 665 else 666 call itrCorrectSclr (y, ac, solinc(1,isclr+5)) 667 endif 668 if (ilset.eq.2 .and. isclr.eq.2) then 669 if (ivconstraint .eq. 1) then 670 call itrBCSclr ( y, ac, iBC, BC, iper, 671 & ilwork) 672c 673c ... applying the volume constraint on second level set scalar 674c 675 call solvecon (y, x, iBC, BC, 676 & iper, ilwork, shp, shgl) 677c 678 endif ! end of volume constraint calculations 679 endif ! end of redistance calculations 680c 681 call itrBCSclr ( y, ac, iBC, BC, iper, 682 & ilwork) 683 endif ! end of flow or scalar update 684 endif ! end of switch between solve or update 685 enddo ! loop over sequence in step 686c 687c 688c.... obtain the time average statistics 689c 690 if (ioform .eq. 2) then 691 692 call stsGetStats( y, yold, ac, acold, 693 & u, uold, x, 694 & shp, shgl, shpb, shglb, 695 & iBC, BC, iper, ilwork, 696 & rowp, colm, lhsK, lhsP ) 697 endif 698 699c 700c Find the solution at the end of the timestep and move it to old 701c 702c 703c ...First to reassign the parameters for the original time integrator scheme 704c 705 if((iLSet.eq.2).and.(ilss.eq.1)) then 706 alfi =alfit 707 gami =gamit 708 almi =almit 709 Delt(1)=Deltt 710 Dtgl =Dtglt 711 endif 712 call itrUpdate( yold, acold, uold, y, ac, u) 713 call itrBC (yold, acold, iBC, BC, iper,ilwork) 714 715 istep = istep + 1 716 lstep = lstep + 1 717c 718c .. Print memory consumption on BGQ 719c 720 call printmeminfo("itrdrv"//char(0)) 721 722c 723c .. Compute vorticity 724c 725 if ( icomputevort == 1) then 726 727 ! vorticity components and magnitude 728 vorticity(:,1) = GradV(:,8)-GradV(:,6) !omega_x 729 vorticity(:,2) = GradV(:,3)-GradV(:,7) !omega_y 730 vorticity(:,3) = GradV(:,4)-GradV(:,2) !omega_z 731 vorticity(:,4) = sqrt( vorticity(:,1)*vorticity(:,1) 732 & + vorticity(:,2)*vorticity(:,2) 733 & + vorticity(:,3)*vorticity(:,3) ) 734 ! Q 735 strain(:,1) = GradV(:,1) !S11 736 strain(:,2) = 0.5*(GradV(:,2)+GradV(:,4)) !S12 737 strain(:,3) = 0.5*(GradV(:,3)+GradV(:,7)) !S13 738 strain(:,4) = GradV(:,5) !S22 739 strain(:,5) = 0.5*(GradV(:,6)+GradV(:,8)) !S23 740 strain(:,6) = GradV(:,9) !S33 741 742 vorticity(:,5) = 0.25*( vorticity(:,4)*vorticity(:,4) !Q 743 & - 2.0*( strain(:,1)*strain(:,1) 744 & + 2* strain(:,2)*strain(:,2) 745 & + 2* strain(:,3)*strain(:,3) 746 & + strain(:,4)*strain(:,4) 747 & + 2* strain(:,5)*strain(:,5) 748 & + strain(:,6)*strain(:,6))) 749 750 endif 751 752c 753c .. write out the solution 754c 755!MR CHANGE 756 757! if (((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) .or. 758! & lstep.eq.nstep(itseq)) then 759 if (((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) .or. 760 & istep.eq.nstep(itseq)) then 761!MR CHANGE END 762 763! Call to restar() will open restart file in write mode (and not append mode) 764! that is needed as other fields are written in append mode 765 call restar ('out ', yold ,ac) 766 if(ideformwall == 1) then 767 call write_displ(myrank, lstep, nshg, 3, uold ) 768 endif 769 770 if(ivort == 1) then 771 call write_field(myrank,'a','vorticity',9,vorticity, 772 & 'd',nshg,5,lstep) 773 endif 774!MR CHANGE 775! if(ioybar.eq.1) then 776! call write_field(myrank,'a','ybar',4, 777! & ybar,'d',nshg,13,lstep) 778! call write_field(myrank,'a','ybar',4, 779! & ybar,'d',nshg,12,lstep) 780! endif 781 782 call printmeminfo("itrdrv after checkpoint"//char(0)) 783!MR CHANGE END 784 785 endif 786 787c 788c ... update the flow history for the impedance convolution, filter it and write it out 789c 790 if(numImpSrfs.gt.zero) then 791 call UpdHistConv(y,nsrflistImp,numImpSrfs) !uses Delt(1) 792 endif 793 794c 795c ... update the flow history for the RCR convolution 796c 797 if(numRCRSrfs.gt.zero) then 798 call UpdHistConv(y,nsrflistRCR,numRCRSrfs) !uses lstep 799 endif 800 801c 802c.... compute the consistent boundary flux 803c 804 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 805 call Bflux ( yold, acold, uold, x, 806 & shp, shgl, shpb, 807 & shglb, ilwork, iBC, 808 & BC, iper, wallssVec) 809 endif 810 811c... dump TIME SERIES 812 813 if (exts) then 814 if (mod(lstep-1,freq).eq.0) then 815 816 if (numpe > 1) then 817 do jj = 1, ntspts 818 vartssoln((jj-1)*ndof+1:jj*ndof)=varts(jj,:) 819 ivarts=zero 820 enddo 821 do k=1,ndof*ntspts 822 if(vartssoln(k).ne.zero) ivarts(k)=1 823 enddo 824 825! call MPI_REDUCE(vartssoln, vartssolng, ndof*ntspts, 826! & MPI_DOUBLE_PRECISION, MPI_SUM, master, 827! & MPI_COMM_WORLD, ierr) 828 829 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 830 call MPI_ALLREDUCE(vartssoln, vartssolng, 831 & ndof*ntspts, 832 & MPI_DOUBLE_PRECISION, MPI_SUM, 833 & MPI_COMM_WORLD, ierr) 834 835! call MPI_REDUCE(ivarts, ivartsg, ndof*ntspts, 836! & MPI_INTEGER, MPI_SUM, master, 837! & MPI_COMM_WORLD, ierr) 838 839 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 840 call MPI_ALLREDUCE(ivarts, ivartsg, ndof*ntspts, 841 & MPI_INTEGER, MPI_SUM, 842 & MPI_COMM_WORLD, ierr) 843 844! if (myrank.eq.zero) then 845 do jj = 1, ntspts 846 847 if(myrank .eq. iv_rank(jj)) then 848 ! No need to update all varts components, only the one treated by the expected rank 849 ! Note: keep varts as a vector, as multiple probes could be treated by one rank 850 indxvarts = (jj-1)*ndof 851 do k=1,ndof 852 if(ivartsg(indxvarts+k).ne.0) then ! none of the vartssoln(parts) were non zero 853 varts(jj,k)=vartssolng(indxvarts+k)/ 854 & ivartsg(indxvarts+k) 855 endif 856 enddo 857 endif !only if myrank eq iv_rank(jj) 858 enddo 859! endif !only on master 860 endif !only if numpe > 1 861 862! if (myrank.eq.zero) then 863 do jj = 1, ntspts 864 if(myrank .eq. iv_rank(jj)) then 865 ifile = 1000+jj 866 write(ifile,555) lstep, (varts(jj,k),k=1,ndof) !Beware of format 555 - check ndof 867c call flush(ifile) 868 if (((irs .ge. 1) .and. 869 & (mod(lstep, ntout) .eq. 0))) then 870 close(ifile) 871 fvarts='varts/varts' 872 fvarts=trim(fvarts)//trim(cname2(jj)) 873 fvarts=trim(fvarts)//trim(cname2(lskeep)) 874 fvarts=trim(fvarts)//'.dat' 875 fvarts=trim(fvarts) 876 open(unit=ifile, file=fvarts, 877 & position='append') 878 endif !only when dumping restart 879 endif 880 enddo 881! endif !only on master 882 883 varts(:,:) = zero ! reset the array for next step 884 885 886!555 format(i6,5(2x,E12.5e2)) 887555 format(i6,6(2x,E20.12e2)) !assuming ndof = 6 here 888 889 endif 890 endif 891 892c 893c.... update and the aerodynamic forces 894c 895 call forces ( yold, ilwork ) 896 897 if((irscale.ge.0).or.(itwmod.gt.0)) 898 & call getvel (yold, ilwork, iBC, 899 & nsons, ifath, velbar) 900 901 if((irscale.ge.0).and.(myrank.eq.master)) then 902 call genscale(yold, x, iper, 903 & iBC, ifath, velbar, 904 & nsons) 905 endif 906c 907c.... print out results. 908c 909 ntoutv=max(ntout,100) ! velb is not needed so often 910 if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then 911 if( (mod(lstep, ntoutv) .eq. 0) .and. 912 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 913 & ((nsonmax.eq.1).and.(iLES.gt.0)))) 914 & call rwvelb ('out ', velbar ,ifail) 915 endif 916c 917c.... end of the NSTEP and NTSEQ loops 918c 919 920 921c 922c.... -------------------> error calculation <----------------- 923c 924 if(ierrcalc.eq.1 .or. ioybar.eq.1) then 925c$$$c 926c$$$c compute average 927c$$$c 928c$$$ tfact=one/istep 929c$$$ ybar =tfact*yold + (one-tfact)*ybar 930 931c compute average 932c ybar(:,1:3) are average velocity components 933c ybar(:,4) is average pressure 934c ybar(:,5) is average speed 935c ybar(:,6:8) is average of sq. of each vel. component 936c ybar(:,9) is average of sq. of pressure 937c ybar(:,10:12) is average of cross vel. components : uv, uw and vw 938c averaging procedure justified only for identical time step sizes 939c ybar(:,13) is average of eddy viscosity 940c ybar(:,14:16) is average vorticity components 941c ybar(:,17) is average vorticity magnitude 942c istep is number of time step 943c 944 icollectybar = 0 945 if(nphasesincycle.eq.0 .or. 946 & istep.gt.ncycles_startphaseavg*nstepsincycle) then 947 icollectybar = 1 948 if((istep-1).eq.ncycles_startphaseavg*nstepsincycle) 949 & istepsinybar = 0 ! init. to zero in first cycle in avg. 950 endif 951 952 if(icollectybar.eq.1) then 953 istepsinybar = istepsinybar+1 954 tfact=one/istepsinybar 955 956 if(myrank.eq.master .and. nphasesincycle.ne.0 .and. 957 & mod((istep-1),nstepsincycle).eq.0) 958 & write(*,*)'nsamples in phase average:',istepsinybar 959 960c ybar to contain the averaged ((u,v,w),p)-fields 961c and speed average, i.e., sqrt(u^2+v^2+w^2) 962c and avg. of sq. terms including 963c u^2, v^2, w^2, p^2 and cross terms of uv, uw and vw 964 965 ybar(:,1) = tfact*yold(:,1) + (one-tfact)*ybar(:,1) 966 ybar(:,2) = tfact*yold(:,2) + (one-tfact)*ybar(:,2) 967 ybar(:,3) = tfact*yold(:,3) + (one-tfact)*ybar(:,3) 968 ybar(:,4) = tfact*yold(:,4) + (one-tfact)*ybar(:,4) 969 ybar(:,5) = tfact*sqrt(yold(:,1)**2+yold(:,2)**2+ 970 & yold(:,3)**2) + (one-tfact)*ybar(:,5) 971 ybar(:,6) = tfact*yold(:,1)**2 + 972 & (one-tfact)*ybar(:,6) 973 ybar(:,7) = tfact*yold(:,2)**2 + 974 & (one-tfact)*ybar(:,7) 975 ybar(:,8) = tfact*yold(:,3)**2 + 976 & (one-tfact)*ybar(:,8) 977 ybar(:,9) = tfact*yold(:,4)**2 + 978 & (one-tfact)*ybar(:,9) 979 ybar(:,10) = tfact*yold(:,1)*yold(:,2) + !uv 980 & (one-tfact)*ybar(:,10) 981 ybar(:,11) = tfact*yold(:,1)*yold(:,3) + !uw 982 & (one-tfact)*ybar(:,11) 983 ybar(:,12) = tfact*yold(:,2)*yold(:,3) + !vw 984 & (one-tfact)*ybar(:,12) 985!MR CHANGE 986 if(nsclr.gt.0) !nut 987 & ybar(:,13) = tfact*yold(:,6) + (one-tfact)*ybar(:,13) 988 989 if(ivort == 1) then !vorticity 990 ybar(:,14) = tfact*vorticity(:,1) + 991 & (one-tfact)*ybar(:,14) 992 ybar(:,15) = tfact*vorticity(:,2) + 993 & (one-tfact)*ybar(:,15) 994 ybar(:,16) = tfact*vorticity(:,3) + 995 & (one-tfact)*ybar(:,16) 996 ybar(:,17) = tfact*vorticity(:,4) + 997 & (one-tfact)*ybar(:,17) 998 endif 999 1000 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1001 wallssVecBar(:,1) = tfact*wallssVec(:,1) 1002 & +(one-tfact)*wallssVecBar(:,1) 1003 wallssVecBar(:,2) = tfact*wallssVec(:,2) 1004 & +(one-tfact)*wallssVecBar(:,2) 1005 wallssVecBar(:,3) = tfact*wallssVec(:,3) 1006 & +(one-tfact)*wallssVecBar(:,3) 1007 endif 1008!MR CHANGE END 1009 endif 1010c 1011c compute phase average 1012c 1013 if(nphasesincycle.ne.0 .and. 1014 & istep.gt.ncycles_startphaseavg*nstepsincycle) then 1015 1016c beginning of cycle is considered as ncycles_startphaseavg*nstepsincycle+1 1017 if((istep-1).eq.ncycles_startphaseavg*nstepsincycle) 1018 & icyclesinavg = 0 ! init. to zero in first cycle in avg. 1019 1020 ! find number of steps between phases 1021 nstepsbtwphase = nstepsincycle/nphasesincycle ! integer value 1022 if(mod(istep-1,nstepsincycle).eq.0) then 1023 iphase = 1 ! init. to one in beginning of every cycle 1024 icyclesinavg = icyclesinavg + 1 1025 endif 1026 1027 icollectphase = 0 1028 istepincycle = mod(istep,nstepsincycle) 1029 if(istepincycle.eq.0) istepincycle=nstepsincycle 1030 if(istepincycle.eq.iphase*nstepsbtwphase) then 1031 icollectphase = 1 1032 iphase = iphase+1 ! use 'iphase-1' below 1033 endif 1034 1035 if(icollectphase.eq.1) then 1036 tfactphase = one/icyclesinavg 1037 1038 if(myrank.eq.master) then 1039 write(*,*) 'nsamples in phase ',iphase-1,': ', 1040 & icyclesinavg 1041 endif 1042 1043 yphbar(:,1,iphase-1) = tfactphase*yold(:,1) + 1044 & (one-tfactphase)*yphbar(:,1,iphase-1) 1045 yphbar(:,2,iphase-1) = tfactphase*yold(:,2) + 1046 & (one-tfactphase)*yphbar(:,2,iphase-1) 1047 yphbar(:,3,iphase-1) = tfactphase*yold(:,3) + 1048 & (one-tfactphase)*yphbar(:,3,iphase-1) 1049 yphbar(:,4,iphase-1) = tfactphase*yold(:,4) + 1050 & (one-tfactphase)*yphbar(:,4,iphase-1) 1051 yphbar(:,5,iphase-1) = tfactphase*sqrt(yold(:,1)**2 1052 & +yold(:,2)**2+yold(:,3)**2) + 1053 & (one-tfactphase)*yphbar(:,5,iphase-1) 1054!MR CHANGE 1055 yphbar(:,6,iphase-1) = 1056 & tfactphase*yold(:,1)*yold(:,1) 1057 & +(one-tfactphase)*yphbar(:,6,iphase-1) 1058 1059 yphbar(:,7,iphase-1) = 1060 & tfactphase*yold(:,1)*yold(:,2) 1061 & +(one-tfactphase)*yphbar(:,7,iphase-1) 1062 1063 yphbar(:,8,iphase-1) = 1064 & tfactphase*yold(:,1)*yold(:,3) 1065 & +(one-tfactphase)*yphbar(:,8,iphase-1) 1066 1067 yphbar(:,9,iphase-1) = 1068 & tfactphase*yold(:,2)*yold(:,2) 1069 & +(one-tfactphase)*yphbar(:,9,iphase-1) 1070 1071 yphbar(:,10,iphase-1) = 1072 & tfactphase*yold(:,2)*yold(:,3) 1073 & +(one-tfactphase)*yphbar(:,10,iphase-1) 1074 1075 yphbar(:,11,iphase-1) = 1076 & tfactphase*yold(:,3)*yold(:,3) 1077 & +(one-tfactphase)*yphbar(:,11,iphase-1) 1078 1079 if(ivort == 1) then 1080 yphbar(:,12,iphase-1) = 1081 & tfactphase*vorticity(:,1) 1082 & +(one-tfactphase)*yphbar(:,12,iphase-1) 1083 yphbar(:,13,iphase-1) = 1084 & tfactphase*vorticity(:,2) 1085 & +(one-tfactphase)*yphbar(:,13,iphase-1) 1086 yphbar(:,14,iphase-1) = 1087 & tfactphase*vorticity(:,3) 1088 & +(one-tfactphase)*yphbar(:,14,iphase-1) 1089 yphbar(:,15,iphase-1) = 1090 & tfactphase*vorticity(:,4) 1091 & +(one-tfactphase)*yphbar(:,15,iphase-1) 1092 endif 1093!MR CHANGE END 1094 endif 1095 endif 1096c 1097c compute rms 1098c 1099 if(icollectybar.eq.1) then 1100 rerr(:, 7)=rerr(:, 7)+(yold(:,1)-ybar(:,1))**2 1101 rerr(:, 8)=rerr(:, 8)+(yold(:,2)-ybar(:,2))**2 1102 rerr(:, 9)=rerr(:, 9)+(yold(:,3)-ybar(:,3))**2 1103 rerr(:,10)=rerr(:,10)+(yold(:,4)-ybar(:,4))**2 1104 endif 1105 endif 1106 1107 if(istop.eq.1000) exit ! stop when delta small (see rstatic) 1108 2000 continue 1109 2001 continue 1110 1111 1112CAD tcorecp2 = second(0) 1113CAD tcorewc2 = second(-1) 1114 1115CAD write(6,*) 'T(core) cpu-wallclock = ',tcorecp2-tcorecp1, 1116CAD & tcorewc2-tcorewc1 1117 1118 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1119 if(myrank.eq.0) then 1120 tcorecp2 = TMRC() 1121 write(6,*) 'T(core) cpu = ',tcorecp2-tcorecp1 1122 write(6,*) '(Elm. form.',tcorecp(1), 1123 & ',Lin. alg. sol.',tcorecp(2),')' 1124 write(6,*) '(Elm. form. Scal.',tcorecpscal(1), 1125 & ',Lin. alg. sol. Scal.',tcorecpscal(2),')' 1126 write(6,*) '' 1127 1128 endif 1129 1130 call print_system_stats(tcorecp, tcorecpscal) 1131 call print_mesh_stats() 1132 call print_mpi_stats() 1133 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1134! return 1135c call MPI_Finalize() 1136c call MPI_ABORT(MPI_COMM_WORLD, ierr) 1137 1138 3000 continue 1139 1140c 1141c.... ----------------------> Post Processing <---------------------- 1142c 1143c.... print out the last step 1144c 1145 if ((irs .ge. 1) .and. ((mod(lstep, ntout) .ne. 0) .or. 1146 & (nstp .eq. 0))) then 1147 if( 1148 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 1149 & ((nsonmax.eq.1).and.iLES.gt.0))) 1150 & call rwvelb ('out ', velbar ,ifail) 1151 endif 1152 1153 1154 lesId = numeqns(1) 1155 1156!MR CHANGE 1157 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1158 if(myrank.eq.0) then 1159 tcormr1 = TMRC() 1160 endif 1161!MR CHANGE END 1162 1163 call saveLesRestart( lesId, aperm , nshg, myrank, lstep, 1164 & nPermDims ) 1165 1166!MR CHANGE 1167 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1168 if(myrank.eq.0) then 1169 tcormr2 = TMRC() 1170 write(6,*) 'Time to call saveLesRestart for projection and'// 1171 & 'pressure projection vectors', tcormr2-tcormr1 1172 endif 1173!MR CHANGE END 1174 1175 1176 if(ierrcalc.eq.1) then 1177c 1178c.....smooth the error indicators 1179c 1180 do i=1,ierrsmooth 1181 call errsmooth( rerr, x, iper, ilwork, shp, shgl, iBC ) 1182 end do 1183 1184!MR CHANGE 1185 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1186 if(myrank.eq.0) then 1187 tcormr1 = TMRC() 1188 endif 1189!MR CHANGE END 1190 1191 call write_error(myrank, lstep, nshg, 10, rerr ) 1192 1193!MR CHANGE 1194 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1195 if(myrank.eq.0) then 1196 tcormr2 = TMRC() 1197 write(6,*) 'Time to write the error fields to the disks', 1198 & tcormr2-tcormr1 1199 endif 1200!MR CHANGE END 1201 1202 1203 endif 1204 1205 if(ioybar.eq.1) then 1206 1207!MR CHANGE 1208! call write_field(myrank,'a','ybar',4, 1209! & ybar,'d',nshg,12,lstep) 1210 if(ivort == 1) then 1211 call write_field(myrank,'a','ybar',4, 1212 & ybar,'d',nshg,17,lstep) 1213 else 1214 call write_field(myrank,'a','ybar',4, 1215 & ybar,'d',nshg,13,lstep) 1216 endif 1217 deallocate(ybar) 1218 1219 1220 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1221 call write_field(myrank,'a','wssbar',6, 1222 & wallssVecBar,'d',nshg,3,lstep) 1223 deallocate(wallssVecbar) 1224 endif 1225 1226!MR CHANGE END 1227 if(nphasesincycle .gt. 0) then 1228 1229 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1230 if(myrank.eq.0) then 1231 tcormr1 = TMRC() 1232 endif 1233 1234 do iphase=1,nphasesincycle 1235 1236! call write_phavg(myrank,'w','phase average',13,iphase, 1237! & yphbar(:,:,iphase),'d',nshg,5,lstep) 1238! ! ybar field is repeated in files for phase average 1239! call write_phavg(myrank,'a','ybar',4,iphase, 1240! & ybar(:,1:5),'d',nshg,5,lstep) 1241 if(ivort == 1) then 1242 call write_phavg2(myrank,'a','phase_average',13,iphase, 1243 & nphasesincycle,yphbar(:,:,iphase),'d',nshg,15,lstep) 1244 else 1245 call write_phavg2(myrank,'a','phase_average',13,iphase, 1246 & nphasesincycle,yphbar(:,:,iphase),'d',nshg,11,lstep) 1247 endif 1248 1249 end do 1250 1251 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1252 if(myrank.eq.0) then 1253 tcormr2 = TMRC() 1254 write(6,*) 'Time to write all phase avg to the disks = ', 1255 & tcormr2-tcormr1 1256 endif 1257 deallocate(yphbar) 1258 endif 1259 1260 endif 1261 1262 if(ivort == 1) then 1263 deallocate(strain,vorticity) 1264 endif 1265 1266 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1267 deallocate(wallssVec) 1268 endif 1269 1270!MR CHANGE END 1271 1272 if ( ( ihessian .eq. 1 ) .and. ( numpe < 2 ) )then 1273 uhess = zero 1274 gradu = zero 1275 tf = zero 1276 1277 do ku=1,nshg 1278c tf(ku,1) = x(ku,1)**2+2*x(ku,1)*x(ku,2) 1279 tf(ku,1) = x(ku,1)**3 1280 end do 1281 1282 call hessian( yold, x, shp, shgl, iBC, 1283 & shpb, shglb, iper, ilwork, uhess, gradu ) 1284 1285 call write_hessian( uhess, gradu, nshg ) 1286 endif 1287 1288c if(iRANS.lt.0 .and. idistcalc.eq.1) then 1289 if(iRANS.lt.0) then 1290c call write_field(myrank,'a','DESd',4, 1291c & elDw,'d',numel,1,lstep) 1292 1293!MR CHANGE 1294 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1295 if(myrank.eq.0) then 1296 tcormr1 = TMRC() 1297 endif 1298!MR CHANGE END 1299 1300 call write_field(myrank,'a','dwal',4,d2wall,'d',nshg,1,lstep) 1301 deallocate(d2wall) 1302 1303!MR CHANGE 1304 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1305 if(myrank.eq.0) then 1306 tcormr2 = TMRC() 1307 write(6,*) 'Time to write dwal to the disks = ', 1308 & tcormr2-tcormr1 1309 endif 1310!MR CHANGE END 1311 1312 1313 endif 1314 1315c 1316c.... close history and aerodynamic forces files 1317c 1318 if (myrank .eq. master) then 1319! close (ihist) 1320 close (iforce) 1321 close(76) 1322 if(numImpSrfs.gt.0 .or. numRCRSrfs.gt.0) then 1323 close (8177) 1324 endif 1325 endif 1326c 1327c.... close varts file for probes 1328c 1329 if(exts) then 1330 do jj=1,ntspts 1331 if (myrank == iv_rank(jj)) then 1332 close(1000+jj) 1333 endif 1334 enddo 1335 deallocate (ivarts) 1336 deallocate (ivartsg) 1337 deallocate (iv_rank) 1338 deallocate (vartssoln) 1339 deallocate (vartssolng) 1340 endif 1341 1342!MR CHANGE 1343 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1344 if(myrank.eq.0) then 1345 write(*,*) 'itrdrv - done with aerodynamic forces' 1346 endif 1347!MR CHANGE 1348 1349 do isrf = 0,MAXSURF 1350! if ( nsrflist(isrf).ne.0 ) then 1351 if ( nsrflist(isrf).ne.0 .and. 1352 & myrank.eq.irankfilesforce(isrf)) then 1353 iunit=60+isrf 1354 close(iunit) 1355 endif 1356 enddo 1357 1358!MR CHANGE 1359 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1360 if(myrank.eq.0) then 1361 write(*,*) 'itrdrv - done with MAXSURF' 1362 endif 1363!MR CHANGE 1364 1365 1366 5 format(1X,F15.10,3X,F15.10,3X,F15.10,3X,F15.10) 1367 444 format(6(2x,e14.7)) 1368c 1369c.... end 1370c 1371 if(nsolflow.eq.1) then 1372 deallocate (lhsK) 1373 deallocate (lhsP) 1374 deallocate (aperm) 1375 deallocate (atemp) 1376 endif 1377 if(nsclrsol.gt.0) then 1378 deallocate (lhsS) 1379 deallocate (apermS) 1380 deallocate (atempS) 1381 endif 1382 1383 if(iabc==1) deallocate(acs) 1384 1385!MR CHANGE 1386 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1387 if(myrank.eq.0) then 1388 write(*,*) 'itrdrv - done - BACK TO process.f' 1389 endif 1390!MR CHANGE 1391 1392 1393 1394 return 1395 end 1396 1397 subroutine lesmodels(y, ac, shgl, shp, 1398 & iper, ilwork, rowp, colm, 1399 & nsons, ifath, x, 1400 & iBC, BC) 1401 1402 include "common.h" 1403 1404 real*8 y(nshg,ndof), ac(nshg,ndof), 1405 & x(numnp,nsd), 1406 & BC(nshg,ndofBC) 1407 real*8 shp(MAXTOP,maxsh,MAXQPT), 1408 & shgl(MAXTOP,nsd,maxsh,MAXQPT) 1409 1410c 1411 integer rowp(nshg,nnz), colm(nshg+1), 1412 & iBC(nshg), 1413 & ilwork(nlwork), 1414 & iper(nshg) 1415 dimension ifath(numnp), nsons(nfath) 1416 1417 real*8, allocatable, dimension(:) :: fwr2,fwr3,fwr4 1418 real*8, allocatable, dimension(:) :: stabdis,cdelsq1 1419 real*8, allocatable, dimension(:,:) :: xavegt, xavegt2,xavegt3 1420 1421 if( (iLES.gt.1) ) then ! Allocate Stuff for advanced LES models 1422 allocate (fwr2(nshg)) 1423 allocate (fwr3(nshg)) 1424 allocate (fwr4(nshg)) 1425 allocate (xavegt(nfath,12)) 1426 allocate (xavegt2(nfath,12)) 1427 allocate (xavegt3(nfath,12)) 1428 allocate (stabdis(nfath)) 1429 endif 1430 1431c.... get dynamic model coefficient 1432c 1433 ilesmod=iLES/10 1434c 1435c digit bit set filter rule, 10 bit set model 1436c 1437 if (ilesmod.eq.0) then ! 0 < iLES< 10 => dyn. model calculated 1438 ! at nodes based on discrete filtering 1439 1440 1441 if(isubmod.eq.2) then 1442 call SUPGdis(y, ac, shgl, shp, 1443 & iper, ilwork, 1444 & nsons, ifath, x, 1445 & iBC, BC, stabdis, xavegt3) 1446 endif 1447 1448 if( ((isubmod.eq.0).or.(isubmod.eq.2)))then ! If no 1449 ! sub-model 1450 ! or SUPG 1451 ! model wanted 1452 1453 if(i2filt.eq.0)then ! If simple filter 1454 1455 if(modlstats .eq. 0) then ! If no model stats wanted 1456 call getdmc (y, shgl, shp, 1457 & iper, ilwork, nsons, 1458 & ifath, x) 1459 else ! else get model stats 1460 call stdfdmc (y, shgl, shp, 1461 & iper, ilwork, nsons, 1462 & ifath, x) 1463 endif ! end of stats if statement 1464 1465 else ! else if twice filtering 1466 1467 call widefdmc(y, shgl, shp, 1468 & iper, ilwork, nsons, 1469 & ifath, x) 1470 1471 1472 endif ! end of simple filter if statement 1473 1474 endif ! end of SUPG or no sub-model if statement 1475 1476 1477 if( (isubmod.eq.1) ) then ! If DFWR sub-model wanted 1478 call cdelBHsq (y, shgl, shp, 1479 & iper, ilwork, nsons, 1480 & ifath, x, cdelsq1) 1481 call FiltRat (y, shgl, shp, 1482 & iper, ilwork, nsons, 1483 & ifath, x, cdelsq1, 1484 & fwr4, fwr3) 1485 1486 1487 if (i2filt.eq.0) then ! If simple filter wanted 1488 call DFWRsfdmc(y, shgl, shp, 1489 & iper, ilwork, nsons, 1490 & ifath, x, fwr2, fwr3) 1491 else ! else if twice filtering wanted 1492 call DFWRwfdmc(y, shgl, shp, 1493 & iper, ilwork, nsons, 1494 & ifath, x, fwr4, fwr4) 1495 endif ! end of simple filter if statement 1496 1497 endif ! end of DFWR sub-model if statement 1498 1499 if( (isubmod.eq.2) )then ! If SUPG sub-model wanted 1500 call dmcSUPG (y, ac, shgl, 1501 & shp, iper, ilwork, 1502 & nsons, ifath, x, 1503 & iBC, BC, rowp, colm, 1504 & xavegt2, stabdis) 1505 endif 1506 1507 if(idis.eq.1)then ! If SUPG/Model dissipation wanted 1508 call ediss (y, ac, shgl, 1509 & shp, iper, ilwork, 1510 & nsons, ifath, x, 1511 & iBC, BC, xavegt) 1512 endif 1513 1514 endif ! end of ilesmod 1515 1516 if (ilesmod .eq. 1) then ! 10 < iLES < 20 => dynamic-mixed 1517 ! at nodes based on discrete filtering 1518 call bardmc (y, shgl, shp, 1519 & iper, ilwork, 1520 & nsons, ifath, x) 1521 endif 1522 1523 if (ilesmod .eq. 2) then ! 20 < iLES < 30 => dynamic at quad 1524 ! pts based on lumped projection filt. 1525 1526 if(isubmod.eq.0)then 1527 call projdmc (y, shgl, shp, 1528 & iper, ilwork, x) 1529 else 1530 call cpjdmcnoi (y, shgl, shp, 1531 & iper, ilwork, x, 1532 & rowp, colm, 1533 & iBC, BC) 1534 endif 1535 1536 endif 1537 1538 if( (iLES.gt.1) ) then ! Deallocate Stuff for advanced LES models 1539 deallocate (fwr2) 1540 deallocate (fwr3) 1541 deallocate (fwr4) 1542 deallocate (xavegt) 1543 deallocate (xavegt2) 1544 deallocate (xavegt3) 1545 deallocate (stabdis) 1546 endif 1547 return 1548 end 1549 1550c 1551c...initialize the coefficients for the impedance convolution 1552c 1553 subroutine CalcImpConvCoef (numISrfs, numTpoints) 1554 1555 use convolImpFlow !uses flow history and impedance for convolution 1556 1557 include "common.h" !for alfi 1558 1559 integer numISrfs, numTpoints 1560 1561 allocate (ConvCoef(numTpoints+2,3)) !same time discret. for all imp. BC 1562 do j=1,numTpoints+2 1563 ConvCoef(j,:)=0.5/numTpoints !dt/2 divided by period T=N*dt 1564 ConvCoef(j,1)=ConvCoef(j,1)*(1.0-alfi)*(1.0-alfi) 1565 ConvCoef(j,2)=ConvCoef(j,2)*(1.0+2*alfi*(1.0-alfi)) 1566 ConvCoef(j,3)=ConvCoef(j,3)*alfi*alfi 1567 enddo 1568 ConvCoef(1,2)=zero 1569 ConvCoef(1,3)=zero 1570 ConvCoef(2,3)=zero 1571 ConvCoef(numTpoints+1,1)=zero 1572 ConvCoef(numTpoints+2,2)=zero 1573 ConvCoef(numTpoints+2,1)=zero 1574c 1575c...calculate the coefficients for the impedance convolution 1576c 1577 allocate (ImpConvCoef(numTpoints+2,numISrfs)) 1578 1579c..coefficients below assume Q linear in time step, Z constant 1580c do j=3,numTpoints 1581c ImpConvCoef(j,:) = ValueListImp(j-1,:)*ConvCoef(j,3) 1582c & + ValueListImp(j,:)*ConvCoef(j,2) 1583c & + ValueListImp(j+1,:)*ConvCoef(j,1) 1584c enddo 1585c ImpConvCoef(1,:) = ValueListImp(2,:)*ConvCoef(1,1) 1586c ImpConvCoef(2,:) = ValueListImp(2,:)*ConvCoef(2,2) 1587c & + ValueListImp(3,:)*ConvCoef(2,1) 1588c ImpConvCoef(numTpoints+1,:) = 1589c & ValueListImp(numTpoints,:)*ConvCoef(numTpoints+1,3) 1590c & + ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+1,2) 1591c ImpConvCoef(numTpoints+2,:) = 1592c & ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+2,3) 1593 1594c..try easiest convolution Q and Z constant per time step 1595 do j=3,numTpoints+1 1596 ImpConvCoef(j,:) = ValueListImp(j-1,:)/numTpoints 1597 enddo 1598 ImpConvCoef(1,:) =zero 1599 ImpConvCoef(2,:) =zero 1600 ImpConvCoef(numTpoints+2,:) = 1601 & ValueListImp(numTpoints+1,:)/numTpoints 1602c compensate for yalpha passed not y in Elmgmr() 1603 ImpConvCoef(numTpoints+1,:)= ImpConvCoef(numTpoints+1,:) 1604 & - ImpConvCoef(numTpoints+2,:)*(1.0-alfi)/alfi 1605 ImpConvCoef(numTpoints+2,:)= ImpConvCoef(numTpoints+2,:)/alfi 1606 return 1607 end 1608 1609c 1610c ... update the flow rate history for the impedance convolution, filter it and write it out 1611c 1612 subroutine UpdHistConv(y,nsrfIdList,numSrfs) 1613 1614 use convolImpFlow !brings ntimeptpT, QHistImp, QHistTry, QHistTryF, numImpSrfs 1615 use convolRCRFlow !brings QHistRCR, numRCRSrfs 1616 1617 include "common.h" 1618 1619 integer nsrfIdList(0:MAXSURF), numSrfs 1620 character*20 fname1 1621 character*10 cname2 1622 character*5 cname 1623 real*8 y(nshg,3) !velocity at time n+1 1624 real*8 NewQ(0:MAXSURF) !temporary unknown for the flow rate 1625 !that needs to be added to the flow history 1626 1627 call GetFlowQ(NewQ,y,nsrfIdList,numSrfs) !new flow at time n+1 1628c 1629c... for imp BC: shift QHist, add new constribution, filter and write out 1630c 1631 if(numImpSrfs.gt.zero) then 1632 do j=1, ntimeptpT 1633 QHistImp(j,1:numSrfs)=QHistImp(j+1,1:numSrfs) 1634 enddo 1635 QHistImp(ntimeptpT+1,1:numSrfs) = NewQ(1:numSrfs) 1636 1637c 1638c....filter the flow rate history 1639c 1640 cutfreq = 10 !hardcoded cutting frequency of the filter 1641 do j=1, ntimeptpT 1642 QHistTry(j,:)=QHistImp(j+1,:) 1643 enddo 1644 call Filter(QHistTryF,QHistTry,ntimeptpT,Delt(1),cutfreq) 1645c.... if no filter applied then uncomment next three lines 1646c do j=1, ntimeptpT 1647c QHistTryF(j,:)=QHistTry(j,:) 1648c enddo 1649 1650c QHistImp(1,:)=zero ! why do we do this? for beta(1,:) = zero it does not really matters 1651 do j=1, ntimeptpT 1652 QHistImp(j+1,:)=QHistTryF(j,:) 1653 enddo 1654c 1655c.... write out the new history of flow rates to Qhistor.dat 1656c 1657 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1658 & (istep .eq. nstep(1)))) .and. 1659 & (myrank .eq. master)) then 1660 open(unit=816, file='Qhistor.dat',status='replace') 1661 write(816,*) ntimeptpT 1662 do j=1,ntimeptpT+1 1663 write(816,*) (QHistImp(j,n),n=1, numSrfs) 1664 enddo 1665 close(816) 1666c... write out a copy with step number to be able to restart 1667 fname1 = 'Qhistor' 1668 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1669 open(unit=8166,file=trim(fname1),status='unknown') 1670 write(8166,*) ntimeptpT 1671 do j=1,ntimeptpT+1 1672 write(8166,*) (QHistImp(j,n),n=1, numSrfs) 1673 enddo 1674 close(8166) 1675 endif 1676 endif 1677 1678c 1679c... for RCR bc just add the new contribution 1680c 1681 if(numRCRSrfs.gt.zero) then 1682 QHistRCR(lstep+1,1:numSrfs) = NewQ(1:numSrfs) 1683c 1684c.... write out the new history of flow rates to Qhistor.dat 1685c 1686 if ((irs .ge. 1) .and. (myrank .eq. master)) then 1687 if(istep.eq.1) then 1688 open(unit=816,file='Qhistor.dat',status='unknown') 1689 else 1690 open(unit=816,file='Qhistor.dat',position='append') 1691 endif 1692 if(istep.eq.1) then 1693 do j=1,lstep 1694 write(816,*) j, (QHistRCR(j,n),n=1,numSrfs) ! read from file of previous run 1695 enddo 1696 endif 1697 write(816,*) lstep+1, (QHistRCR(lstep+1,n),n=1, numSrfs) 1698 close(816) 1699c... write out a copy with step number to be able to restart 1700 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1701 & (istep .eq. nstep(1)))) .and. 1702 & (myrank .eq. master)) then 1703 fname1 = 'Qhistor' 1704 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1705 open(unit=8166,file=trim(fname1),status='unknown') 1706 write(8166,*) lstep+1 1707 do j=1,lstep+1 1708 write(8166,*) (QHistRCR(j,n),n=1, numSrfs) 1709 enddo 1710 close(8166) 1711 endif 1712 endif 1713 endif 1714 1715 return 1716 end 1717 1718c 1719c...calculate the time varying coefficients for the RCR convolution 1720c 1721 subroutine CalcRCRConvCoef (stepn, numSrfs) 1722 1723 use convolRCRFlow !brings in ValueListRCR, dtRCR 1724 1725 include "common.h" !brings alfi 1726 1727 integer numSrfs, stepn 1728 1729 RCRConvCoef = zero 1730 if (stepn .eq. 0) then 1731 RCRConvCoef(1,:) = ValueListRCR(1,:)*(1.0-alfi) + 1732 & ValueListRCR(3,:)*(-alfi + 1.0 + 1/dtRCR(:) 1733 & - exp(-alfi*dtRCR(:))*(1 + 1/dtRCR(:))) 1734 RCRConvCoef(2,:) = ValueListRCR(1,:)*alfi 1735 & + ValueListRCR(3,:) 1736 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1737 endif 1738 if (stepn .ge. 1) then 1739 RCRConvCoef(1,:) =-ValueListRCR(3,:)*exp(-dtRCR(:)*(stepn+alfi)) 1740 & *(1 + (1 - exp(dtRCR(:)))/dtRCR(:)) 1741 RCRConvCoef(stepn+1,:) = ValueListRCR(1,:)*(1-alfi) 1742 & - ValueListRCR(3,:)*(alfi - 1 - 1/dtRCR(:) 1743 & + exp(-alfi*dtRCR(:))/dtRCR(:)*(2 - exp(-dtRCR(:)))) 1744 RCRConvCoef(stepn+2,:) = ValueListRCR(1,:)*alfi 1745 & + ValueListRCR(3,:) 1746 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1747 endif 1748 if (stepn .ge. 2) then 1749 do j=2,stepn 1750 RCRConvCoef(j,:) = ValueListRCR(3,:)/dtRCR(:)* 1751 & exp(-dtRCR(:)*(stepn + alfi + 2 - j))* 1752 & (1 - exp(dtRCR(:)))**2 1753 enddo 1754 endif 1755 1756c compensate for yalpha passed not y in Elmgmr() 1757 RCRConvCoef(stepn+1,:)= RCRConvCoef(stepn+1,:) 1758 & - RCRConvCoef(stepn+2,:)*(1.0-alfi)/alfi 1759 RCRConvCoef(stepn+2,:)= RCRConvCoef(stepn+2,:)/alfi 1760 1761 return 1762 end 1763 1764c 1765c...calculate the time dependent H operator for the RCR convolution 1766c 1767 subroutine CalcHopRCR (timestepRCR, stepn, numSrfs) 1768 1769 use convolRCRFlow !brings in HopRCR, dtRCR 1770 1771 include "common.h" 1772 1773 integer numSrfs, stepn 1774 real*8 PdistCur(0:MAXSURF), timestepRCR 1775 1776 HopRCR=zero 1777 call RCRint(timestepRCR*(stepn + alfi),PdistCur) 1778 HopRCR(1:numSrfs) = RCRic(1:numSrfs) 1779 & *exp(-dtRCR(1:numSrfs)*(stepn + alfi)) + PdistCur(1:numSrfs) 1780 return 1781 end 1782c 1783c ... initialize the influence of the initial conditions for the RCR BC 1784c 1785 subroutine calcRCRic(y,srfIdList,numSrfs) 1786 1787 use convolRCRFlow !brings RCRic, ValueListRCR, ValuePdist 1788 1789 include "common.h" 1790 1791 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled 1792 real*8 y(nshg,4) !need velocity and pressure 1793 real*8 Qini(0:MAXSURF) !initial flow rate 1794 real*8 PdistIni(0:MAXSURF) !initial distal pressure 1795 real*8 Pini(0:MAXSURF),CoupleArea(0:MAXSURF) ! initial pressure 1796 real*8 VelOnly(nshg,3), POnly(nshg) 1797 1798 allocate (RCRic(0:MAXSURF)) 1799 1800 if(lstep.eq.0) then 1801 VelOnly(:,1:3)=y(:,1:3) 1802 call GetFlowQ(Qini,VelOnly,srfIdList,numSrfs) !get initial flow 1803 QHistRCR(1,1:numSrfs)=Qini(1:numSrfs) !initialize QHistRCR 1804 POnly(:)=y(:,4) ! pressure 1805 call integrScalar(Pini,POnly,srfIdList,numSrfs) !get initial pressure integral 1806 POnly(:)=one ! one to get area 1807 call integrScalar(CoupleArea,POnly,srfIdList,numSrfs) !get surf area 1808 Pini(1:numSrfs) = Pini(1:numSrfs)/CoupleArea(1:numSrfs) 1809 else 1810 Qini(1:numSrfs)=QHistRCR(1,1:numSrfs) 1811 Pini(1:numSrfs)=zero ! hack 1812 endif 1813 call RCRint(istep,PdistIni) !get initial distal P (use istep) 1814 RCRic(1:numSrfs) = Pini(1:numSrfs) 1815 & - ValueListRCR(1,:)*Qini(1:numSrfs)-PdistIni(1:numSrfs) 1816 return 1817 end 1818 1819c.........function that integrates a scalar over a boundary 1820 subroutine integrScalar(scalInt,scal,srfIdList,numSrfs) 1821 1822 use pvsQbi !brings ndsurf, NASC 1823 1824 include "common.h" 1825 include "mpif.h" 1826 1827 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled, i, k 1828 real*8 scal(nshg), scalInt(0:MAXSURF), scalIntProc(0:MAXSURF) 1829 1830 scalIntProc = zero 1831 do i = 1,nshg 1832 if(numSrfs.gt.zero) then 1833 do k = 1,numSrfs 1834 irankCoupled = 0 1835 if (srfIdList(k).eq.ndsurf(i)) then 1836 irankCoupled=k 1837 scalIntProc(irankCoupled) = scalIntProc(irankCoupled) 1838 & + NASC(i)*scal(i) 1839 exit 1840 endif 1841 enddo 1842 endif 1843 enddo 1844c 1845c at this point, each scalint has its "nodes" contributions to the scalar 1846c accumulated into scalIntProc. Note, because NASC is on processor this 1847c will NOT be the scalar for the surface yet 1848c 1849c.... reduce integrated scalar for each surface, push on scalInt 1850c 1851 npars=MAXSURF+1 1852 call MPI_ALLREDUCE (scalIntProc, scalInt(:), npars, 1853 & MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr) 1854 1855 return 1856 end 1857 1858 1859 1860 1861 1862