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