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