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