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