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