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 do 2000 istp = 1, nstp 468 if(iramp.eq.1) 469 & call BCprofileScale(vbc_prof,BC,yold) 470 471 call rerun_check(stopjob) 472 if(myrank.eq.master) write(*,*) 473 & 'stopjob,lstep,istep', stopjob,lstep,istep 474 if(stopjob.eq.lstep) then 475 stopjob=-2 ! this is the code to finish 476 if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then 477 if(myrank.eq.master) write(*,*) 478 & 'line 473 says last step written so exit' 479 goto 2002 ! the step was written last step so just exit 480 else 481 if(myrank.eq.master) 482 & write(*,*) 'line 473 says last step not written' 483 istep=nstp ! have to do this so that solution will be written 484 goto 2001 485 endif 486 endif 487 488 xi=istp*1.0/nstp 489 datmat(1,2,1)=rmub*(1.0-xi)+xi*rmue 490c write(*,*) "current mol. visc = ", datmat(1,2,1) 491c.... if we have time varying boundary conditions update the values of BC. 492c these will be for time step n+1 so use lstep+1 493c 494 if(itvn.gt.0) call BCint((lstep+1)*Delt(1), shp, shgl, 495 & shpb, shglb, x, BC, iBC) 496 497c 498c ... calculate the pressure contribution that depends on the history for the Imp. BC 499c 500 if(numImpSrfs.gt.0) then 501 call pHist(poldImp,QHistImp,ImpConvCoef, 502 & ntimeptpT,numImpSrfs) 503 if (myrank.eq.master) 504 & write(8177,*) (poldImp(n), n=1,numImpSrfs) 505 endif 506c 507c ... calc the pressure contribution that depends on the history for the RCR BC 508c 509 if(numRCRSrfs.gt.0) then 510 call CalcHopRCR (Delt(itseq), lstep, numRCRSrfs) 511 call CalcRCRConvCoef(lstep,numRCRSrfs) 512 call pHist(poldRCR,QHistRCR,RCRConvCoef,nsteprcr, 513 & numRCRSrfs) 514 if (myrank.eq.master) 515 & write(8177,*) (poldRCR(n), n=1,numRCRSrfs) 516 endif 517c 518c Decay of scalars 519c 520 if(nsclr.gt.0 .and. tdecay.ne.1) then 521 yold(:,6:ndof)=y(:,6:ndof)*tdecay 522 BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*tdecay 523 endif 524 525 if(nosource.eq.1) BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*0.8 526 527 528 if(iLES.gt.0) then !complicated stuff has moved to 529 !routine below 530 call lesmodels(yold, acold, shgl, shp, 531 & iper, ilwork, rowp, colm, 532 & nsons, ifath, x, 533 & iBC, BC) 534 535 536 endif 537 538c.... set traction BCs for modeled walls 539c 540 if (itwmod.ne.0) then 541 call asbwmod(yold, acold, x, BC, iBC, 542 & iper, ilwork, ifath, velbar) 543 endif 544 545!MR CHANGE 546c 547c.... Determine whether the vorticity field needs to be computed for this time step or not 548c 549 icomputevort = 0 550 if (ivort == 1) then ! Print vorticity = True in solver.inp 551 ! We then compute the vorticity only if we 552 ! 1) we write an intermediate checkpoint 553 ! 2) we reach the last time step and write the last checkpoint 554 ! 3) we accumulate statistics in ybar for every time step 555 ! BEWARE: we need here lstep+1 and istep+1 because the lstep and 556 ! istep gets incremened after the flowsolve, further below 557 if (((irs .ge. 1) .and. (mod(lstep+1, ntout) .eq. 0)) .or. 558 & istep+1.eq.nstep(itseq) .or. ioybar == 1) then 559 icomputevort = 1 560 endif 561 endif 562 563! write(*,*) 'icomputevort: ',icomputevort, ' - istep: ', 564! & istep,' - nstep(itseq):',nstep(itseq),'- lstep:', 565! & lstep, '- ntout:', ntout 566!MR CHANGE 567 568c 569c.... -----------------------> predictor phase <----------------------- 570c 571 call itrPredict(yold, y, acold, ac , uold, u, iBC) 572 call itrBC (y, ac, iBC, BC, iper,ilwork) 573 574 if(nsolt.eq.1) then 575 isclr=0 576 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 577 endif 578 do isclr=1,nsclr 579 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 580 enddo 581 iter=0 582 ilss=0 ! this is a switch thrown on first solve of LS redistance 583 do istepc=1,seqsize 584 icode=stepseq(istepc) 585 if(mod(icode,5).eq.0) then ! this is a solve 586 isolve=icode/10 587 if(icode.eq.0) then ! flow solve (encoded as 0) 588c 589 iter = iter+1 590 ifuncs(1) = ifuncs(1) + 1 591c 592 Force(1) = zero 593 Force(2) = zero 594 Force(3) = zero 595 HFlux = zero 596 lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 597 598 call SolFlow(y, ac, u, 599 & yold, acold, uold, 600 & x, iBC, 601 & BC, res, 602 & nPermDims, nTmpDims, aperm, 603 & atemp, iper, 604 & ilwork, shp, shgl, 605 & shpb, shglb, rowp, 606 & colm, lhsK, lhsP, 607 & solinc, rerr, tcorecp, 608 & GradV) 609 610 else ! scalar type solve 611 if (icode.eq.5) then ! Solve for Temperature 612 ! (encoded as (nsclr+1)*10) 613 isclr=0 614 ifuncs(2) = ifuncs(2) + 1 615 j=1 616 else ! solve a scalar (encoded at isclr*10) 617 isclr=isolve 618 ifuncs(isclr+2) = ifuncs(isclr+2) + 1 619 j=isclr+nsolt 620 if((iLSet.eq.2).and.(ilss.eq.0) 621 & .and.(isclr.eq.2)) then 622 ilss=1 ! throw switch (once per step) 623 y(:,7)=y(:,6) ! redistance field initialized 624 ac(:,7) = zero 625 call itrBCSclr ( y, ac, iBC, BC, iper, 626 & ilwork) 627c 628c....store the flow alpha, gamma parameter values and assigm them the 629c....Backward Euler parameters to solve the second levelset scalar 630c 631 alfit=alfi 632 gamit=gami 633 almit=almi 634 Deltt=Delt(1) 635 Dtglt=Dtgl 636 alfi = 1 637 gami = 1 638 almi = 1 639c Delt(1)= Deltt ! Give a pseudo time step 640 Dtgl = one / Delt(1) 641 endif ! level set eq. 2 642 endif ! deciding between temperature and scalar 643 644 lhs = 1 - min(1,mod(ifuncs(isclr+2)-1, 645 & LHSupd(isclr+2))) 646 647 call SolSclr(y, ac, u, 648 & yold, acold, uold, 649 & x, iBC, 650 & BC, nPermDimsS,nTmpDimsS, 651 & apermS(1,1,j), atempS, iper, 652 & ilwork, shp, shgl, 653 & shpb, shglb, rowp, 654 & colm, lhsS(1,j), 655 & solinc(1,isclr+5), tcorecpscal) 656 657 658 endif ! end of scalar type solve 659 660 else ! this is an update (mod did not equal zero) 661 iupdate=icode/10 ! what to update 662 if(icode.eq.1) then !update flow 663 call itrCorrect ( y, ac, u, solinc, iBC) 664 call itrBC (y, ac, iBC, BC, iper, ilwork) 665 else ! update scalar 666 isclr=iupdate !unless 667 if(icode.eq.6) isclr=0 668 if(iRANS.lt.-100)then ! RANS 669 call itrCorrectSclrPos(y,ac,solinc(1,isclr+5)) 670 else 671 call itrCorrectSclr (y, ac, solinc(1,isclr+5)) 672 endif 673 if (ilset.eq.2 .and. isclr.eq.2) then 674 if (ivconstraint .eq. 1) then 675 call itrBCSclr ( y, ac, iBC, BC, iper, 676 & ilwork) 677c 678c ... applying the volume constraint on second level set scalar 679c 680 call solvecon (y, x, iBC, BC, 681 & iper, ilwork, shp, shgl) 682c 683 endif ! end of volume constraint calculations 684 endif ! end of redistance calculations 685c 686 call itrBCSclr ( y, ac, iBC, BC, iper, 687 & ilwork) 688 endif ! end of flow or scalar update 689 endif ! end of switch between solve or update 690 enddo ! loop over sequence in step 691c 692c 693c.... obtain the time average statistics 694c 695 if (ioform .eq. 2) then 696 697 call stsGetStats( y, yold, ac, acold, 698 & u, uold, x, 699 & shp, shgl, shpb, shglb, 700 & iBC, BC, iper, ilwork, 701 & rowp, colm, lhsK, lhsP ) 702 endif 703 704c 705c Find the solution at the end of the timestep and move it to old 706c 707c 708c ...First to reassign the parameters for the original time integrator scheme 709c 710 if((iLSet.eq.2).and.(ilss.eq.1)) then 711 alfi =alfit 712 gami =gamit 713 almi =almit 714 Delt(1)=Deltt 715 Dtgl =Dtglt 716 endif 717 call itrUpdate( yold, acold, uold, y, ac, u) 718 call itrBC (yold, acold, iBC, BC, iper,ilwork) 719 720 istep = istep + 1 721 lstep = lstep + 1 722c 723c .. Print memory consumption on BGQ 724c 725 call printmeminfo("itrdrv"//char(0)) 726 727c 728c .. Compute vorticity 729c 730 if ( icomputevort == 1) then 731 732 ! vorticity components and magnitude 733 vorticity(:,1) = GradV(:,8)-GradV(:,6) !omega_x 734 vorticity(:,2) = GradV(:,3)-GradV(:,7) !omega_y 735 vorticity(:,3) = GradV(:,4)-GradV(:,2) !omega_z 736 vorticity(:,4) = sqrt( vorticity(:,1)*vorticity(:,1) 737 & + vorticity(:,2)*vorticity(:,2) 738 & + vorticity(:,3)*vorticity(:,3) ) 739 ! Q 740 strain(:,1) = GradV(:,1) !S11 741 strain(:,2) = 0.5*(GradV(:,2)+GradV(:,4)) !S12 742 strain(:,3) = 0.5*(GradV(:,3)+GradV(:,7)) !S13 743 strain(:,4) = GradV(:,5) !S22 744 strain(:,5) = 0.5*(GradV(:,6)+GradV(:,8)) !S23 745 strain(:,6) = GradV(:,9) !S33 746 747 vorticity(:,5) = 0.25*( vorticity(:,4)*vorticity(:,4) !Q 748 & - 2.0*( strain(:,1)*strain(:,1) 749 & + 2* strain(:,2)*strain(:,2) 750 & + 2* strain(:,3)*strain(:,3) 751 & + strain(:,4)*strain(:,4) 752 & + 2* strain(:,5)*strain(:,5) 753 & + strain(:,6)*strain(:,6))) 754 755 endif 756c 757c .. write out the instantaneoud solution 758c 7592001 continue ! we could get here by 2001 label if user requested stop 760 if (((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) .or. 761 & istep.eq.nstep(itseq)) then 762 763!so that we can see progress in force file close it so that it flushes 764!and then reopen in append mode 765 766 close(iforce) 767 open (unit=iforce, file=fforce, position='append') 768 769! Call to restar() will open restart file in write mode (and not append mode) 770! that is needed as other fields are written in append mode 771 772 call restar ('out ', yold ,ac) 773 if(ideformwall == 1) then 774 call write_displ(myrank, lstep, nshg, 3, uold ) 775 endif 776 777 if(ivort == 1) then 778 call write_field(myrank,'a','vorticity',9,vorticity, 779 & 'd',nshg,5,lstep) 780 endif 781 782 call printmeminfo("itrdrv after checkpoint"//char(0)) 783 else if(stopjob.eq.-2) then 784 if(myrank.eq.master) then 785 write(*,*) 'line 755 says no write before stopping' 786 write(*,*) 'istep,nstep,irs',istep,nstep(itseq),irs 787 endif 788 endif !just the instantaneous stuff for vidos 789c 790c.... compute the consistent boundary flux 791c 792 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 793 call Bflux ( yold, acold, uold, x, 794 & shp, shgl, shpb, 795 & shglb, ilwork, iBC, 796 & BC, iper, wallssVec) 797 endif 798 799 if(stopjob.eq.-2) goto 2003 800 801 802c 803c ... update the flow history for the impedance convolution, filter it and write it out 804c 805 if(numImpSrfs.gt.zero) then 806 call UpdHistConv(y,nsrflistImp,numImpSrfs) !uses Delt(1) 807 endif 808 809c 810c ... update the flow history for the RCR convolution 811c 812 if(numRCRSrfs.gt.zero) then 813 call UpdHistConv(y,nsrflistRCR,numRCRSrfs) !uses lstep 814 endif 815 816 817c... dump TIME SERIES 818 819 if (exts) then 820 if (mod(lstep-1,freq).eq.0) then 821 822 if (numpe > 1) then 823 do jj = 1, ntspts 824 vartssoln((jj-1)*ndof+1:jj*ndof)=varts(jj,:) 825 ivarts=zero 826 enddo 827 do k=1,ndof*ntspts 828 if(vartssoln(k).ne.zero) ivarts(k)=1 829 enddo 830 831! call MPI_REDUCE(vartssoln, vartssolng, ndof*ntspts, 832! & MPI_DOUBLE_PRECISION, MPI_SUM, master, 833! & MPI_COMM_WORLD, ierr) 834 835 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 836 call MPI_ALLREDUCE(vartssoln, vartssolng, 837 & ndof*ntspts, 838 & MPI_DOUBLE_PRECISION, MPI_SUM, 839 & MPI_COMM_WORLD, ierr) 840 841! call MPI_REDUCE(ivarts, ivartsg, ndof*ntspts, 842! & MPI_INTEGER, MPI_SUM, master, 843! & MPI_COMM_WORLD, ierr) 844 845 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 846 call MPI_ALLREDUCE(ivarts, ivartsg, ndof*ntspts, 847 & MPI_INTEGER, MPI_SUM, 848 & MPI_COMM_WORLD, ierr) 849 850! if (myrank.eq.zero) then 851 do jj = 1, ntspts 852 853 if(myrank .eq. iv_rank(jj)) then 854 ! No need to update all varts components, only the one treated by the expected rank 855 ! Note: keep varts as a vector, as multiple probes could be treated by one rank 856 indxvarts = (jj-1)*ndof 857 do k=1,ndof 858 if(ivartsg(indxvarts+k).ne.0) then ! none of the vartssoln(parts) were non zero 859 varts(jj,k)=vartssolng(indxvarts+k)/ 860 & ivartsg(indxvarts+k) 861 endif 862 enddo 863 endif !only if myrank eq iv_rank(jj) 864 enddo 865! endif !only on master 866 endif !only if numpe > 1 867 868! if (myrank.eq.zero) then 869 do jj = 1, ntspts 870 if(myrank .eq. iv_rank(jj)) then 871 ifile = 1000+jj 872 write(ifile,555) lstep, (varts(jj,k),k=1,ndof) !Beware of format 555 - check ndof 873c call flush(ifile) 874 if (((irs .ge. 1) .and. 875 & (mod(lstep, ntout) .eq. 0))) then 876 close(ifile) 877 fvarts='varts/varts' 878 fvarts=trim(fvarts)//trim(cname2(jj)) 879 fvarts=trim(fvarts)//trim(cname2(lskeep)) 880 fvarts=trim(fvarts)//'.dat' 881 fvarts=trim(fvarts) 882 open(unit=ifile, file=fvarts, 883 & position='append') 884 endif !only when dumping restart 885 endif 886 enddo 887! endif !only on master 888 889 varts(:,:) = zero ! reset the array for next step 890 891 892!555 format(i6,5(2x,E12.5e2)) 893555 format(i6,6(2x,E20.12e2)) !assuming ndof = 6 here 894 895 endif 896 endif 897 898c 899c.... update and the aerodynamic forces 900c 901 call forces ( yold, ilwork ) 902 903 if((irscale.ge.0).or.(itwmod.gt.0)) 904 & call getvel (yold, ilwork, iBC, 905 & nsons, ifath, velbar) 906 907 if((irscale.ge.0).and.(myrank.eq.master)) then 908 call genscale(yold, x, iper, 909 & iBC, ifath, velbar, 910 & nsons) 911 endif 912c 913c.... print out results. 914c 915 ntoutv=max(ntout,100) ! velb is not needed so often 916 if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then 917 if( (mod(lstep, ntoutv) .eq. 0) .and. 918 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 919 & ((nsonmax.eq.1).and.(iLES.gt.0)))) 920 & call rwvelb ('out ', velbar ,ifail) 921 endif 922c 923c.... end of the NSTEP and NTSEQ loops 924c 925 926 927c 928c.... -------------------> error calculation <----------------- 929c 930 if(ierrcalc.eq.1 .or. ioybar.eq.1) then 931c$$$c 932c$$$c compute average 933c$$$c 934c$$$ tfact=one/istep 935c$$$ ybar =tfact*yold + (one-tfact)*ybar 936 937c compute average 938c ybar(:,1:3) are average velocity components 939c ybar(:,4) is average pressure 940c ybar(:,5) is average speed 941c ybar(:,6:8) is average of sq. of each vel. component 942c ybar(:,9) is average of sq. of pressure 943c ybar(:,10:12) is average of cross vel. components : uv, uw and vw 944c averaging procedure justified only for identical time step sizes 945c ybar(:,13) is average of eddy viscosity 946c ybar(:,14:16) is average vorticity components 947c ybar(:,17) is average vorticity magnitude 948c istep is number of time step 949c 950 icollectybar = 0 951 if(nphasesincycle.eq.0 .or. 952 & istep.gt.ncycles_startphaseavg*nstepsincycle) then 953 icollectybar = 1 954 if((istep-1).eq.ncycles_startphaseavg*nstepsincycle) 955 & istepsinybar = 0 ! init. to zero in first cycle in avg. 956 endif 957 958 if(icollectybar.eq.1) then 959 istepsinybar = istepsinybar+1 960 tfact=one/istepsinybar 961 962 if(myrank.eq.master .and. nphasesincycle.ne.0 .and. 963 & mod((istep-1),nstepsincycle).eq.0) 964 & write(*,*)'nsamples in phase average:',istepsinybar 965 966c ybar to contain the averaged ((u,v,w),p)-fields 967c and speed average, i.e., sqrt(u^2+v^2+w^2) 968c and avg. of sq. terms including 969c u^2, v^2, w^2, p^2 and cross terms of uv, uw and vw 970 971 ybar(:,1) = tfact*yold(:,1) + (one-tfact)*ybar(:,1) 972 ybar(:,2) = tfact*yold(:,2) + (one-tfact)*ybar(:,2) 973 ybar(:,3) = tfact*yold(:,3) + (one-tfact)*ybar(:,3) 974 ybar(:,4) = tfact*yold(:,4) + (one-tfact)*ybar(:,4) 975 ybar(:,5) = tfact*sqrt(yold(:,1)**2+yold(:,2)**2+ 976 & yold(:,3)**2) + (one-tfact)*ybar(:,5) 977 ybar(:,6) = tfact*yold(:,1)**2 + 978 & (one-tfact)*ybar(:,6) 979 ybar(:,7) = tfact*yold(:,2)**2 + 980 & (one-tfact)*ybar(:,7) 981 ybar(:,8) = tfact*yold(:,3)**2 + 982 & (one-tfact)*ybar(:,8) 983 ybar(:,9) = tfact*yold(:,4)**2 + 984 & (one-tfact)*ybar(:,9) 985 ybar(:,10) = tfact*yold(:,1)*yold(:,2) + !uv 986 & (one-tfact)*ybar(:,10) 987 ybar(:,11) = tfact*yold(:,1)*yold(:,3) + !uw 988 & (one-tfact)*ybar(:,11) 989 ybar(:,12) = tfact*yold(:,2)*yold(:,3) + !vw 990 & (one-tfact)*ybar(:,12) 991!MR CHANGE 992 if(nsclr.gt.0) !nut 993 & ybar(:,13) = tfact*yold(:,6) + (one-tfact)*ybar(:,13) 994 995 if(ivort == 1) then !vorticity 996 ybar(:,14) = tfact*vorticity(:,1) + 997 & (one-tfact)*ybar(:,14) 998 ybar(:,15) = tfact*vorticity(:,2) + 999 & (one-tfact)*ybar(:,15) 1000 ybar(:,16) = tfact*vorticity(:,3) + 1001 & (one-tfact)*ybar(:,16) 1002 ybar(:,17) = tfact*vorticity(:,4) + 1003 & (one-tfact)*ybar(:,17) 1004 endif 1005 1006 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1007 wallssVecBar(:,1) = tfact*wallssVec(:,1) 1008 & +(one-tfact)*wallssVecBar(:,1) 1009 wallssVecBar(:,2) = tfact*wallssVec(:,2) 1010 & +(one-tfact)*wallssVecBar(:,2) 1011 wallssVecBar(:,3) = tfact*wallssVec(:,3) 1012 & +(one-tfact)*wallssVecBar(:,3) 1013 endif 1014!MR CHANGE END 1015 endif 1016c 1017c compute phase average 1018c 1019 if(nphasesincycle.ne.0 .and. 1020 & istep.gt.ncycles_startphaseavg*nstepsincycle) then 1021 1022c beginning of cycle is considered as ncycles_startphaseavg*nstepsincycle+1 1023 if((istep-1).eq.ncycles_startphaseavg*nstepsincycle) 1024 & icyclesinavg = 0 ! init. to zero in first cycle in avg. 1025 1026 ! find number of steps between phases 1027 nstepsbtwphase = nstepsincycle/nphasesincycle ! integer value 1028 if(mod(istep-1,nstepsincycle).eq.0) then 1029 iphase = 1 ! init. to one in beginning of every cycle 1030 icyclesinavg = icyclesinavg + 1 1031 endif 1032 1033 icollectphase = 0 1034 istepincycle = mod(istep,nstepsincycle) 1035 if(istepincycle.eq.0) istepincycle=nstepsincycle 1036 if(istepincycle.eq.iphase*nstepsbtwphase) then 1037 icollectphase = 1 1038 iphase = iphase+1 ! use 'iphase-1' below 1039 endif 1040 1041 if(icollectphase.eq.1) then 1042 tfactphase = one/icyclesinavg 1043 1044 if(myrank.eq.master) then 1045 write(*,*) 'nsamples in phase ',iphase-1,': ', 1046 & icyclesinavg 1047 endif 1048 1049 yphbar(:,1,iphase-1) = tfactphase*yold(:,1) + 1050 & (one-tfactphase)*yphbar(:,1,iphase-1) 1051 yphbar(:,2,iphase-1) = tfactphase*yold(:,2) + 1052 & (one-tfactphase)*yphbar(:,2,iphase-1) 1053 yphbar(:,3,iphase-1) = tfactphase*yold(:,3) + 1054 & (one-tfactphase)*yphbar(:,3,iphase-1) 1055 yphbar(:,4,iphase-1) = tfactphase*yold(:,4) + 1056 & (one-tfactphase)*yphbar(:,4,iphase-1) 1057 yphbar(:,5,iphase-1) = tfactphase*sqrt(yold(:,1)**2 1058 & +yold(:,2)**2+yold(:,3)**2) + 1059 & (one-tfactphase)*yphbar(:,5,iphase-1) 1060!MR CHANGE 1061 yphbar(:,6,iphase-1) = 1062 & tfactphase*yold(:,1)*yold(:,1) 1063 & +(one-tfactphase)*yphbar(:,6,iphase-1) 1064 1065 yphbar(:,7,iphase-1) = 1066 & tfactphase*yold(:,1)*yold(:,2) 1067 & +(one-tfactphase)*yphbar(:,7,iphase-1) 1068 1069 yphbar(:,8,iphase-1) = 1070 & tfactphase*yold(:,1)*yold(:,3) 1071 & +(one-tfactphase)*yphbar(:,8,iphase-1) 1072 1073 yphbar(:,9,iphase-1) = 1074 & tfactphase*yold(:,2)*yold(:,2) 1075 & +(one-tfactphase)*yphbar(:,9,iphase-1) 1076 1077 yphbar(:,10,iphase-1) = 1078 & tfactphase*yold(:,2)*yold(:,3) 1079 & +(one-tfactphase)*yphbar(:,10,iphase-1) 1080 1081 yphbar(:,11,iphase-1) = 1082 & tfactphase*yold(:,3)*yold(:,3) 1083 & +(one-tfactphase)*yphbar(:,11,iphase-1) 1084 1085 if(ivort == 1) then 1086 yphbar(:,12,iphase-1) = 1087 & tfactphase*vorticity(:,1) 1088 & +(one-tfactphase)*yphbar(:,12,iphase-1) 1089 yphbar(:,13,iphase-1) = 1090 & tfactphase*vorticity(:,2) 1091 & +(one-tfactphase)*yphbar(:,13,iphase-1) 1092 yphbar(:,14,iphase-1) = 1093 & tfactphase*vorticity(:,3) 1094 & +(one-tfactphase)*yphbar(:,14,iphase-1) 1095 yphbar(:,15,iphase-1) = 1096 & tfactphase*vorticity(:,4) 1097 & +(one-tfactphase)*yphbar(:,15,iphase-1) 1098 endif 1099!MR CHANGE END 1100 endif 1101 endif 1102c 1103c compute rms 1104c 1105 if(icollectybar.eq.1) then 1106 rerr(:, 7)=rerr(:, 7)+(yold(:,1)-ybar(:,1))**2 1107 rerr(:, 8)=rerr(:, 8)+(yold(:,2)-ybar(:,2))**2 1108 rerr(:, 9)=rerr(:, 9)+(yold(:,3)-ybar(:,3))**2 1109 rerr(:,10)=rerr(:,10)+(yold(:,4)-ybar(:,4))**2 1110 endif 1111 endif 1112 2003 continue ! we get here if stopjob equals lstep and this jumped over 1113! the statistics computation because we have no new data to average in 1114! rather we are just trying to output the last state that was not already 1115! written 1116c 1117c.... ----------------------> Complete Restart Processing <---------------------- 1118c 1119! for now it is the same frequency but need to change this 1120! soon.... but don't forget to change the field counter in 1121! new_interface.cc 1122! 1123 if (((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) .or. 1124 & istep.eq.nstep(itseq)) then 1125 if ((irs .ge. 1) .and. ((mod(lstep, ntout) .ne. 0) .or. 1126 & (nstp .eq. 0))) then 1127 if( 1128 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 1129 & ((nsonmax.eq.1).and.iLES.gt.0))) 1130 & call rwvelb ('out ', velbar ,ifail) 1131 endif 1132 1133 lesId = numeqns(1) 1134 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1135 if(myrank.eq.0) then 1136 tcormr1 = TMRC() 1137 endif 1138 call saveLesRestart( lesId, aperm , nshg, myrank, lstep, 1139 & nPermDims ) 1140 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1141 if(myrank.eq.0) then 1142 tcormr2 = TMRC() 1143 write(6,*) 'call saveLesRestart for projection and'// 1144 & 'pressure projection vectors', tcormr2-tcormr1 1145 endif 1146 1147 if(ierrcalc.eq.1) then 1148c 1149c.....smooth the error indicators 1150c 1151 do i=1,ierrsmooth 1152 call errsmooth( rerr, x, iper, ilwork, shp, shgl, iBC ) 1153 end do 1154 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1155 if(myrank.eq.0) then 1156 tcormr1 = TMRC() 1157 endif 1158 call write_error(myrank, lstep, nshg, 10, rerr ) 1159 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1160 if(myrank.eq.0) then 1161 tcormr2 = TMRC() 1162 write(6,*) 'Time to write the error fields to the disks', 1163 & tcormr2-tcormr1 1164 endif 1165 endif ! ierrcalc 1166 1167 if(ioybar.eq.1) then 1168 if(ivort == 1) then 1169 call write_field(myrank,'a','ybar',4, 1170 & ybar,'d',nshg,17,lstep) 1171 else 1172 call write_field(myrank,'a','ybar',4, 1173 & ybar,'d',nshg,13,lstep) 1174 endif 1175 1176 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1177 call write_field(myrank,'a','wssbar',6, 1178 & wallssVecBar,'d',nshg,3,lstep) 1179 endif 1180 1181 if(nphasesincycle .gt. 0) then 1182 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1183 if(myrank.eq.0) then 1184 tcormr1 = TMRC() 1185 endif 1186 do iphase=1,nphasesincycle 1187 if(ivort == 1) then 1188 call write_phavg2(myrank,'a','phase_average',13,iphase, 1189 & nphasesincycle,yphbar(:,:,iphase),'d',nshg,15,lstep) 1190 else 1191 call write_phavg2(myrank,'a','phase_average',13,iphase, 1192 & nphasesincycle,yphbar(:,:,iphase),'d',nshg,11,lstep) 1193 endif 1194 end do 1195 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1196 if(myrank.eq.0) then 1197 tcormr2 = TMRC() 1198 write(6,*) 'write all phase avg to the disks = ', 1199 & tcormr2-tcormr1 1200 endif 1201 endif !nphasesincyle 1202 endif !ioybar 1203 1204 if ( ( ihessian .eq. 1 ) .and. ( numpe < 2 ) )then 1205 uhess = zero 1206 gradu = zero 1207 tf = zero 1208 1209 do ku=1,nshg 1210 tf(ku,1) = x(ku,1)**3 1211 end do 1212 call hessian( yold, x, shp, shgl, iBC, 1213 & shpb, shglb, iper, ilwork, uhess, gradu ) 1214 1215 call write_hessian( uhess, gradu, nshg ) 1216 endif 1217 1218 if(iRANS.lt.0) then 1219 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1220 if(myrank.eq.0) then 1221 tcormr1 = TMRC() 1222 endif 1223 call write_field(myrank,'a','dwal',4,d2wall,'d', 1224 & nshg,1,lstep) 1225 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1226 if(myrank.eq.0) then 1227 tcormr2 = TMRC() 1228 write(6,*) 'Time to write dwal to the disks = ', 1229 & tcormr2-tcormr1 1230 endif 1231 endif !iRANS 1232 1233 endif ! write out complete restart state 1234 !next 2 lines are two ways to end early 1235 if(stopjob.eq.-2) goto 2002 1236 if(istop.eq.1000) goto 2002 ! stop when delta small (see rstatic) 1237 2000 continue 1238 2002 continue 1239 1240! dnoe with time stepping so deallocate fields alfready written 1241! 1242 if(ioybar.eq.1) then 1243 deallocate(ybar) 1244 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1245 deallocate(wallssVecbar) 1246 endif 1247 if(nphasesincycle .gt. 0) then 1248 deallocate(yphbar) 1249 endif !nphasesincyle 1250 endif !ioybar 1251 if(ivort == 1) then 1252 deallocate(strain,vorticity) 1253 endif 1254 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1255 deallocate(wallssVec) 1256 endif 1257 if(iRANS.lt.0) then 1258 deallocate(d2wall) 1259 endif 1260 1261 1262CAD tcorecp2 = second(0) 1263CAD tcorewc2 = second(-1) 1264 1265CAD write(6,*) 'T(core) cpu-wallclock = ',tcorecp2-tcorecp1, 1266CAD & tcorewc2-tcorewc1 1267 1268 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1269 if(myrank.eq.0) then 1270 tcorecp2 = TMRC() 1271 write(6,*) 'T(core) cpu = ',tcorecp2-tcorecp1 1272 write(6,*) '(Elm. form.',tcorecp(1), 1273 & ',Lin. alg. sol.',tcorecp(2),')' 1274 write(6,*) '(Elm. form. Scal.',tcorecpscal(1), 1275 & ',Lin. alg. sol. Scal.',tcorecpscal(2),')' 1276 write(6,*) '' 1277 1278 endif 1279 1280 call print_system_stats(tcorecp, tcorecpscal) 1281 call print_mesh_stats() 1282 call print_mpi_stats() 1283 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1284! return 1285c call MPI_Finalize() 1286c call MPI_ABORT(MPI_COMM_WORLD, ierr) 1287 1288 3000 continue 1289 1290 1291c 1292c.... close history and aerodynamic forces files 1293c 1294 if (myrank .eq. master) then 1295! close (ihist) 1296 close (iforce) 1297 close(76) 1298 if(numImpSrfs.gt.0 .or. numRCRSrfs.gt.0) then 1299 close (8177) 1300 endif 1301 endif 1302c 1303c.... close varts file for probes 1304c 1305 if(exts) then 1306 do jj=1,ntspts 1307 if (myrank == iv_rank(jj)) then 1308 close(1000+jj) 1309 endif 1310 enddo 1311 deallocate (ivarts) 1312 deallocate (ivartsg) 1313 deallocate (iv_rank) 1314 deallocate (vartssoln) 1315 deallocate (vartssolng) 1316 endif 1317 1318!MR CHANGE 1319 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1320 if(myrank.eq.0) then 1321 write(*,*) 'itrdrv - done with aerodynamic forces' 1322 endif 1323!MR CHANGE 1324 1325 do isrf = 0,MAXSURF 1326! if ( nsrflist(isrf).ne.0 ) then 1327 if ( nsrflist(isrf).ne.0 .and. 1328 & myrank.eq.irankfilesforce(isrf)) then 1329 iunit=60+isrf 1330 close(iunit) 1331 endif 1332 enddo 1333 1334!MR CHANGE 1335 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1336 if(myrank.eq.0) then 1337 write(*,*) 'itrdrv - done with MAXSURF' 1338 endif 1339!MR CHANGE 1340 1341 1342 5 format(1X,F15.10,3X,F15.10,3X,F15.10,3X,F15.10) 1343 444 format(6(2x,e14.7)) 1344c 1345c.... end 1346c 1347 if(nsolflow.eq.1) then 1348 deallocate (lhsK) 1349 deallocate (lhsP) 1350 deallocate (aperm) 1351 deallocate (atemp) 1352 endif 1353 if(nsclrsol.gt.0) then 1354 deallocate (lhsS) 1355 deallocate (apermS) 1356 deallocate (atempS) 1357 endif 1358 1359 if(iabc==1) deallocate(acs) 1360 1361!MR CHANGE 1362 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1363 if(myrank.eq.0) then 1364 write(*,*) 'itrdrv - done - BACK TO process.f' 1365 endif 1366!MR CHANGE 1367 1368 1369 1370 return 1371 end 1372 1373 subroutine lesmodels(y, ac, shgl, shp, 1374 & iper, ilwork, rowp, colm, 1375 & nsons, ifath, x, 1376 & iBC, BC) 1377 1378 include "common.h" 1379 1380 real*8 y(nshg,ndof), ac(nshg,ndof), 1381 & x(numnp,nsd), 1382 & BC(nshg,ndofBC) 1383 real*8 shp(MAXTOP,maxsh,MAXQPT), 1384 & shgl(MAXTOP,nsd,maxsh,MAXQPT) 1385 1386c 1387 integer rowp(nshg,nnz), colm(nshg+1), 1388 & iBC(nshg), 1389 & ilwork(nlwork), 1390 & iper(nshg) 1391 dimension ifath(numnp), nsons(nfath) 1392 1393 real*8, allocatable, dimension(:) :: fwr2,fwr3,fwr4 1394 real*8, allocatable, dimension(:) :: stabdis,cdelsq1 1395 real*8, allocatable, dimension(:,:) :: xavegt, xavegt2,xavegt3 1396 1397 if( (iLES.gt.1) ) then ! Allocate Stuff for advanced LES models 1398 allocate (fwr2(nshg)) 1399 allocate (fwr3(nshg)) 1400 allocate (fwr4(nshg)) 1401 allocate (xavegt(nfath,12)) 1402 allocate (xavegt2(nfath,12)) 1403 allocate (xavegt3(nfath,12)) 1404 allocate (stabdis(nfath)) 1405 endif 1406 1407c.... get dynamic model coefficient 1408c 1409 ilesmod=iLES/10 1410c 1411c digit bit set filter rule, 10 bit set model 1412c 1413 if (ilesmod.eq.0) then ! 0 < iLES< 10 => dyn. model calculated 1414 ! at nodes based on discrete filtering 1415 1416 1417 if(isubmod.eq.2) then 1418 call SUPGdis(y, ac, shgl, shp, 1419 & iper, ilwork, 1420 & nsons, ifath, x, 1421 & iBC, BC, stabdis, xavegt3) 1422 endif 1423 1424 if( ((isubmod.eq.0).or.(isubmod.eq.2)))then ! If no 1425 ! sub-model 1426 ! or SUPG 1427 ! model wanted 1428 1429 if(i2filt.eq.0)then ! If simple filter 1430 1431 if(modlstats .eq. 0) then ! If no model stats wanted 1432 call getdmc (y, shgl, shp, 1433 & iper, ilwork, nsons, 1434 & ifath, x) 1435 else ! else get model stats 1436 call stdfdmc (y, shgl, shp, 1437 & iper, ilwork, nsons, 1438 & ifath, x) 1439 endif ! end of stats if statement 1440 1441 else ! else if twice filtering 1442 1443 call widefdmc(y, shgl, shp, 1444 & iper, ilwork, nsons, 1445 & ifath, x) 1446 1447 1448 endif ! end of simple filter if statement 1449 1450 endif ! end of SUPG or no sub-model if statement 1451 1452 1453 if( (isubmod.eq.1) ) then ! If DFWR sub-model wanted 1454 call cdelBHsq (y, shgl, shp, 1455 & iper, ilwork, nsons, 1456 & ifath, x, cdelsq1) 1457 call FiltRat (y, shgl, shp, 1458 & iper, ilwork, nsons, 1459 & ifath, x, cdelsq1, 1460 & fwr4, fwr3) 1461 1462 1463 if (i2filt.eq.0) then ! If simple filter wanted 1464 call DFWRsfdmc(y, shgl, shp, 1465 & iper, ilwork, nsons, 1466 & ifath, x, fwr2, fwr3) 1467 else ! else if twice filtering wanted 1468 call DFWRwfdmc(y, shgl, shp, 1469 & iper, ilwork, nsons, 1470 & ifath, x, fwr4, fwr4) 1471 endif ! end of simple filter if statement 1472 1473 endif ! end of DFWR sub-model if statement 1474 1475 if( (isubmod.eq.2) )then ! If SUPG sub-model wanted 1476 call dmcSUPG (y, ac, shgl, 1477 & shp, iper, ilwork, 1478 & nsons, ifath, x, 1479 & iBC, BC, rowp, colm, 1480 & xavegt2, stabdis) 1481 endif 1482 1483 if(idis.eq.1)then ! If SUPG/Model dissipation wanted 1484 call ediss (y, ac, shgl, 1485 & shp, iper, ilwork, 1486 & nsons, ifath, x, 1487 & iBC, BC, xavegt) 1488 endif 1489 1490 endif ! end of ilesmod 1491 1492 if (ilesmod .eq. 1) then ! 10 < iLES < 20 => dynamic-mixed 1493 ! at nodes based on discrete filtering 1494 call bardmc (y, shgl, shp, 1495 & iper, ilwork, 1496 & nsons, ifath, x) 1497 endif 1498 1499 if (ilesmod .eq. 2) then ! 20 < iLES < 30 => dynamic at quad 1500 ! pts based on lumped projection filt. 1501 1502 if(isubmod.eq.0)then 1503 call projdmc (y, shgl, shp, 1504 & iper, ilwork, x) 1505 else 1506 call cpjdmcnoi (y, shgl, shp, 1507 & iper, ilwork, x, 1508 & rowp, colm, 1509 & iBC, BC) 1510 endif 1511 1512 endif 1513 1514 if( (iLES.gt.1) ) then ! Deallocate Stuff for advanced LES models 1515 deallocate (fwr2) 1516 deallocate (fwr3) 1517 deallocate (fwr4) 1518 deallocate (xavegt) 1519 deallocate (xavegt2) 1520 deallocate (xavegt3) 1521 deallocate (stabdis) 1522 endif 1523 return 1524 end 1525 1526c 1527c...initialize the coefficients for the impedance convolution 1528c 1529 subroutine CalcImpConvCoef (numISrfs, numTpoints) 1530 1531 use convolImpFlow !uses flow history and impedance for convolution 1532 1533 include "common.h" !for alfi 1534 1535 integer numISrfs, numTpoints 1536 1537 allocate (ConvCoef(numTpoints+2,3)) !same time discret. for all imp. BC 1538 do j=1,numTpoints+2 1539 ConvCoef(j,:)=0.5/numTpoints !dt/2 divided by period T=N*dt 1540 ConvCoef(j,1)=ConvCoef(j,1)*(1.0-alfi)*(1.0-alfi) 1541 ConvCoef(j,2)=ConvCoef(j,2)*(1.0+2*alfi*(1.0-alfi)) 1542 ConvCoef(j,3)=ConvCoef(j,3)*alfi*alfi 1543 enddo 1544 ConvCoef(1,2)=zero 1545 ConvCoef(1,3)=zero 1546 ConvCoef(2,3)=zero 1547 ConvCoef(numTpoints+1,1)=zero 1548 ConvCoef(numTpoints+2,2)=zero 1549 ConvCoef(numTpoints+2,1)=zero 1550c 1551c...calculate the coefficients for the impedance convolution 1552c 1553 allocate (ImpConvCoef(numTpoints+2,numISrfs)) 1554 1555c..coefficients below assume Q linear in time step, Z constant 1556c do j=3,numTpoints 1557c ImpConvCoef(j,:) = ValueListImp(j-1,:)*ConvCoef(j,3) 1558c & + ValueListImp(j,:)*ConvCoef(j,2) 1559c & + ValueListImp(j+1,:)*ConvCoef(j,1) 1560c enddo 1561c ImpConvCoef(1,:) = ValueListImp(2,:)*ConvCoef(1,1) 1562c ImpConvCoef(2,:) = ValueListImp(2,:)*ConvCoef(2,2) 1563c & + ValueListImp(3,:)*ConvCoef(2,1) 1564c ImpConvCoef(numTpoints+1,:) = 1565c & ValueListImp(numTpoints,:)*ConvCoef(numTpoints+1,3) 1566c & + ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+1,2) 1567c ImpConvCoef(numTpoints+2,:) = 1568c & ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+2,3) 1569 1570c..try easiest convolution Q and Z constant per time step 1571 do j=3,numTpoints+1 1572 ImpConvCoef(j,:) = ValueListImp(j-1,:)/numTpoints 1573 enddo 1574 ImpConvCoef(1,:) =zero 1575 ImpConvCoef(2,:) =zero 1576 ImpConvCoef(numTpoints+2,:) = 1577 & ValueListImp(numTpoints+1,:)/numTpoints 1578c compensate for yalpha passed not y in Elmgmr() 1579 ImpConvCoef(numTpoints+1,:)= ImpConvCoef(numTpoints+1,:) 1580 & - ImpConvCoef(numTpoints+2,:)*(1.0-alfi)/alfi 1581 ImpConvCoef(numTpoints+2,:)= ImpConvCoef(numTpoints+2,:)/alfi 1582 return 1583 end 1584 1585c 1586c ... update the flow rate history for the impedance convolution, filter it and write it out 1587c 1588 subroutine UpdHistConv(y,nsrfIdList,numSrfs) 1589 1590 use convolImpFlow !brings ntimeptpT, QHistImp, QHistTry, QHistTryF, numImpSrfs 1591 use convolRCRFlow !brings QHistRCR, numRCRSrfs 1592 1593 include "common.h" 1594 1595 integer nsrfIdList(0:MAXSURF), numSrfs 1596 character*20 fname1 1597 character*10 cname2 1598 character*5 cname 1599 real*8 y(nshg,3) !velocity at time n+1 1600 real*8 NewQ(0:MAXSURF) !temporary unknown for the flow rate 1601 !that needs to be added to the flow history 1602 1603 call GetFlowQ(NewQ,y,nsrfIdList,numSrfs) !new flow at time n+1 1604c 1605c... for imp BC: shift QHist, add new constribution, filter and write out 1606c 1607 if(numImpSrfs.gt.zero) then 1608 do j=1, ntimeptpT 1609 QHistImp(j,1:numSrfs)=QHistImp(j+1,1:numSrfs) 1610 enddo 1611 QHistImp(ntimeptpT+1,1:numSrfs) = NewQ(1:numSrfs) 1612 1613c 1614c....filter the flow rate history 1615c 1616 cutfreq = 10 !hardcoded cutting frequency of the filter 1617 do j=1, ntimeptpT 1618 QHistTry(j,:)=QHistImp(j+1,:) 1619 enddo 1620 call Filter(QHistTryF,QHistTry,ntimeptpT,Delt(1),cutfreq) 1621c.... if no filter applied then uncomment next three lines 1622c do j=1, ntimeptpT 1623c QHistTryF(j,:)=QHistTry(j,:) 1624c enddo 1625 1626c QHistImp(1,:)=zero ! why do we do this? for beta(1,:) = zero it does not really matters 1627 do j=1, ntimeptpT 1628 QHistImp(j+1,:)=QHistTryF(j,:) 1629 enddo 1630c 1631c.... write out the new history of flow rates to Qhistor.dat 1632c 1633 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1634 & (istep .eq. nstep(1)))) .and. 1635 & (myrank .eq. master)) then 1636 open(unit=816, file='Qhistor.dat',status='replace') 1637 write(816,*) ntimeptpT 1638 do j=1,ntimeptpT+1 1639 write(816,*) (QHistImp(j,n),n=1, numSrfs) 1640 enddo 1641 close(816) 1642c... write out a copy with step number to be able to restart 1643 fname1 = 'Qhistor' 1644 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1645 open(unit=8166,file=trim(fname1),status='unknown') 1646 write(8166,*) ntimeptpT 1647 do j=1,ntimeptpT+1 1648 write(8166,*) (QHistImp(j,n),n=1, numSrfs) 1649 enddo 1650 close(8166) 1651 endif 1652 endif 1653 1654c 1655c... for RCR bc just add the new contribution 1656c 1657 if(numRCRSrfs.gt.zero) then 1658 QHistRCR(lstep+1,1:numSrfs) = NewQ(1:numSrfs) 1659c 1660c.... write out the new history of flow rates to Qhistor.dat 1661c 1662 if ((irs .ge. 1) .and. (myrank .eq. master)) then 1663 if(istep.eq.1) then 1664 open(unit=816,file='Qhistor.dat',status='unknown') 1665 else 1666 open(unit=816,file='Qhistor.dat',position='append') 1667 endif 1668 if(istep.eq.1) then 1669 do j=1,lstep 1670 write(816,*) j, (QHistRCR(j,n),n=1,numSrfs) ! read from file of previous run 1671 enddo 1672 endif 1673 write(816,*) lstep+1, (QHistRCR(lstep+1,n),n=1, numSrfs) 1674 close(816) 1675c... write out a copy with step number to be able to restart 1676 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1677 & (istep .eq. nstep(1)))) .and. 1678 & (myrank .eq. master)) then 1679 fname1 = 'Qhistor' 1680 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1681 open(unit=8166,file=trim(fname1),status='unknown') 1682 write(8166,*) lstep+1 1683 do j=1,lstep+1 1684 write(8166,*) (QHistRCR(j,n),n=1, numSrfs) 1685 enddo 1686 close(8166) 1687 endif 1688 endif 1689 endif 1690 1691 return 1692 end 1693 1694c 1695c...calculate the time varying coefficients for the RCR convolution 1696c 1697 subroutine CalcRCRConvCoef (stepn, numSrfs) 1698 1699 use convolRCRFlow !brings in ValueListRCR, dtRCR 1700 1701 include "common.h" !brings alfi 1702 1703 integer numSrfs, stepn 1704 1705 RCRConvCoef = zero 1706 if (stepn .eq. 0) then 1707 RCRConvCoef(1,:) = ValueListRCR(1,:)*(1.0-alfi) + 1708 & ValueListRCR(3,:)*(-alfi + 1.0 + 1/dtRCR(:) 1709 & - exp(-alfi*dtRCR(:))*(1 + 1/dtRCR(:))) 1710 RCRConvCoef(2,:) = ValueListRCR(1,:)*alfi 1711 & + ValueListRCR(3,:) 1712 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1713 endif 1714 if (stepn .ge. 1) then 1715 RCRConvCoef(1,:) =-ValueListRCR(3,:)*exp(-dtRCR(:)*(stepn+alfi)) 1716 & *(1 + (1 - exp(dtRCR(:)))/dtRCR(:)) 1717 RCRConvCoef(stepn+1,:) = ValueListRCR(1,:)*(1-alfi) 1718 & - ValueListRCR(3,:)*(alfi - 1 - 1/dtRCR(:) 1719 & + exp(-alfi*dtRCR(:))/dtRCR(:)*(2 - exp(-dtRCR(:)))) 1720 RCRConvCoef(stepn+2,:) = ValueListRCR(1,:)*alfi 1721 & + ValueListRCR(3,:) 1722 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1723 endif 1724 if (stepn .ge. 2) then 1725 do j=2,stepn 1726 RCRConvCoef(j,:) = ValueListRCR(3,:)/dtRCR(:)* 1727 & exp(-dtRCR(:)*(stepn + alfi + 2 - j))* 1728 & (1 - exp(dtRCR(:)))**2 1729 enddo 1730 endif 1731 1732c compensate for yalpha passed not y in Elmgmr() 1733 RCRConvCoef(stepn+1,:)= RCRConvCoef(stepn+1,:) 1734 & - RCRConvCoef(stepn+2,:)*(1.0-alfi)/alfi 1735 RCRConvCoef(stepn+2,:)= RCRConvCoef(stepn+2,:)/alfi 1736 1737 return 1738 end 1739 1740c 1741c...calculate the time dependent H operator for the RCR convolution 1742c 1743 subroutine CalcHopRCR (timestepRCR, stepn, numSrfs) 1744 1745 use convolRCRFlow !brings in HopRCR, dtRCR 1746 1747 include "common.h" 1748 1749 integer numSrfs, stepn 1750 real*8 PdistCur(0:MAXSURF), timestepRCR 1751 1752 HopRCR=zero 1753 call RCRint(timestepRCR*(stepn + alfi),PdistCur) 1754 HopRCR(1:numSrfs) = RCRic(1:numSrfs) 1755 & *exp(-dtRCR(1:numSrfs)*(stepn + alfi)) + PdistCur(1:numSrfs) 1756 return 1757 end 1758c 1759c ... initialize the influence of the initial conditions for the RCR BC 1760c 1761 subroutine calcRCRic(y,srfIdList,numSrfs) 1762 1763 use convolRCRFlow !brings RCRic, ValueListRCR, ValuePdist 1764 1765 include "common.h" 1766 1767 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled 1768 real*8 y(nshg,4) !need velocity and pressure 1769 real*8 Qini(0:MAXSURF) !initial flow rate 1770 real*8 PdistIni(0:MAXSURF) !initial distal pressure 1771 real*8 Pini(0:MAXSURF),CoupleArea(0:MAXSURF) ! initial pressure 1772 real*8 VelOnly(nshg,3), POnly(nshg) 1773 1774 allocate (RCRic(0:MAXSURF)) 1775 1776 if(lstep.eq.0) then 1777 VelOnly(:,1:3)=y(:,1:3) 1778 call GetFlowQ(Qini,VelOnly,srfIdList,numSrfs) !get initial flow 1779 QHistRCR(1,1:numSrfs)=Qini(1:numSrfs) !initialize QHistRCR 1780 POnly(:)=y(:,4) ! pressure 1781 call integrScalar(Pini,POnly,srfIdList,numSrfs) !get initial pressure integral 1782 POnly(:)=one ! one to get area 1783 call integrScalar(CoupleArea,POnly,srfIdList,numSrfs) !get surf area 1784 Pini(1:numSrfs) = Pini(1:numSrfs)/CoupleArea(1:numSrfs) 1785 else 1786 Qini(1:numSrfs)=QHistRCR(1,1:numSrfs) 1787 Pini(1:numSrfs)=zero ! hack 1788 endif 1789 call RCRint(istep,PdistIni) !get initial distal P (use istep) 1790 RCRic(1:numSrfs) = Pini(1:numSrfs) 1791 & - ValueListRCR(1,:)*Qini(1:numSrfs)-PdistIni(1:numSrfs) 1792 return 1793 end 1794 1795c.........function that integrates a scalar over a boundary 1796 subroutine integrScalar(scalInt,scal,srfIdList,numSrfs) 1797 1798 use pvsQbi !brings ndsurf, NASC 1799 1800 include "common.h" 1801 include "mpif.h" 1802 1803 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled, i, k 1804 real*8 scal(nshg), scalInt(0:MAXSURF), scalIntProc(0:MAXSURF) 1805 1806 scalIntProc = zero 1807 do i = 1,nshg 1808 if(numSrfs.gt.zero) then 1809 do k = 1,numSrfs 1810 irankCoupled = 0 1811 if (srfIdList(k).eq.ndsurf(i)) then 1812 irankCoupled=k 1813 scalIntProc(irankCoupled) = scalIntProc(irankCoupled) 1814 & + NASC(i)*scal(i) 1815 exit 1816 endif 1817 enddo 1818 endif 1819 enddo 1820c 1821c at this point, each scalint has its "nodes" contributions to the scalar 1822c accumulated into scalIntProc. Note, because NASC is on processor this 1823c will NOT be the scalar for the surface yet 1824c 1825c.... reduce integrated scalar for each surface, push on scalInt 1826c 1827 npars=MAXSURF+1 1828 call MPI_ALLREDUCE (scalIntProc, scalInt(:), npars, 1829 & MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr) 1830 1831 return 1832 end 1833 1834 1835 1836 1837 1838