1 subroutine itrdrv (y, ac, uold, x, 2 & iBC, BC, 3 & iper, ilwork, shp, 4 & shgl, shpb, shglb, 5 & ifath, velbar, nsons ) 6c 7c---------------------------------------------------------------------- 8c 9c This iterative driver is the semi-discrete, predictor multi-corrector 10c algorithm. It contains the Hulbert Generalized Alpha method which 11c is 2nd order accurate for Rho_inf from 0 to 1. The method can be 12c made first-order accurate by setting Rho_inf=-1. It uses a 13c GMRES iterative solver. 14c 15c working arrays: 16c y (nshg,ndof) : Y variables 17c x (nshg,nsd) : node coordinates 18c iBC (nshg) : BC codes 19c BC (nshg,ndofBC) : BC constraint parameters 20c iper (nshg) : periodicity table 21c 22c shape functions: 23c shp (nshape,ngauss) : interior element shape functions 24c shgl (nsd,nshape,ngauss) : local shape function gradients 25c shpb (nshapeb,ngaussb) : boundary element shape functions 26c shglb (nsd,nshapeb,ngaussb) : bdry. elt. shape gradients 27c 28c Zdenek Johan, Winter 1991. (Fortran 90) 29c---------------------------------------------------------------------- 30c 31 use pvsQbi !gives us splag (the spmass at the end of this run 32 use specialBC !gives us itvn 33 use timedataC !allows collection of time series 34 use MachControl !PID to control the inlet velocity. 35 use blowerControl !gives us BC_enable 36 use turbSA 37 use wallData 38 use fncorpmod 39 40 include "common.h" 41 include "mpif.h" 42 include "auxmpi.h" 43 44c 45 dimension y(nshg,ndof), ac(nshg,ndof), 46 & yold(nshg,ndof), acold(nshg,ndof), 47 & x(numnp,nsd), iBC(nshg), 48 & BC(nshg,ndofBC), ilwork(nlwork), 49 & iper(nshg), uold(nshg,nsd) 50c 51 dimension res(nshg,nflow), 52 & rest(nshg), solinc(nshg,ndof) 53c 54 dimension shp(MAXTOP,maxsh,MAXQPT), 55 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 56 & shpb(MAXTOP,maxsh,MAXQPT), 57 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 58 real*8 almit, alfit, gamit 59 dimension ifath(numnp), velbar(nfath,ndof), nsons(nfath) 60 real*8 rerr(nshg,10),ybar(nshg,ndof+8) ! 8 is for avg. of square as uu, vv, ww, pp, TT, uv, uw, and vw 61 real*8, allocatable, dimension(:,:) :: vortG 62 real*8, allocatable, dimension(:,:,:) :: BDiag 63 64! integer, allocatable, dimension(:) :: iv_rank !used with MRasquin's version of probe points 65! integer :: iv_rankpernode, iv_totnodes, iv_totcores 66! integer :: iv_node, iv_core, iv_thread 67 68! assuming three profiles to control (inlet, bottom FC and top FC) 69! store velocity profile set via BC 70 real*8 vbc_prof(nshg,3) 71 character(len=60) fvarts 72 integer ifuncs(6), iarray(10) 73 real*8 elDw(numel) ! element average of DES d variable 74c 75c Here are the data structures for sparse matrix GMRES 76c 77 integer, allocatable, dimension(:,:) :: rowp 78 integer, allocatable, dimension(:) :: colm 79 real*8, allocatable, dimension(:,:) :: lhsK 80 real*8, allocatable, dimension(:,:) :: EGmass 81 real*8, allocatable, dimension(:,:) :: EGmasst 82 83 integer iTurbWall(nshg) 84 real*8 yInlet(3), yInletg(3) 85 integer imapped, imapInlet(nshg) !for now, used for setting Blower conditions 86! real*8 M_th, M_tc, M_tt 87! logical exMc 88! real*8 vBC, vBCg 89 real*8 vortmax, vortmaxg 90 91 iprec=0 !PETSc - Disable PHASTA's BDiag. TODO: Preprocssor Switch 92 93 call findTurbWall(iTurbWall) 94 95!------- 96! SETUP 97!------- 98 99 !HACK for debugging suction 100! call Write_Debug(myrank, 'wallNormal'//char(0), 101! & 'wnorm'//char(0), wnorm, 102! & 'd', nshg, 3, lstep) 103 104 !Probe Point Setup 105 call initProbePoints() 106 if(exts) then !exts is set in initProbePoints 107 write(fvarts, "('./varts/varts.', I0, '.dat')") lstep 108 fvarts = trim(fvarts) 109 110 if(myrank .eq. master) then 111 call TD_writeHeader(fvarts) 112 endif 113 endif 114 115 !Mach Control Setup 116 call MC_init(Delt, lstep, BC) 117 exMC = exMC .and. exts !If probe points aren't available, turn 118 !the Mach Control off 119 if(exMC) then 120 call MC_applyBC(BC) 121 call MC_printState() 122 endif 123 124 125c 126c.... open history and aerodynamic forces files 127c 128 if (myrank .eq. master) then 129 open (unit=ihist, file=fhist, status='unknown') 130 open (unit=iforce, file=fforce, status='unknown') 131 endif 132c 133c 134c.... initialize 135c 136 ifuncs = 0 ! func. evaluation counter 137 istep = 0 138 ntotGM = 0 ! number of GMRES iterations 139 time = 0 140 yold = y 141 acold = ac 142 143!Blower Setup 144 call BC_init(Delt, lstep, BC) !Note: sets BC_enable 145! fix the yold values to the reset BC 146 if(BC_enable) call itrBC (yold, ac, iBC, BC, iper, ilwork) 147! without the above, second solve of first steps is fouled up 148! 149 if (mod(impl(1),100)/10 .eq. 1) then 150c 151c generate the sparse data fill vectors 152c 153 allocate (rowp(nshg,nnz)) 154 allocate (colm(nshg+1)) 155 call genadj(colm, rowp, icnt ) ! preprocess the adjacency list 156 157 nnz_tot=icnt ! this is exactly the number of non-zero 158 ! blocks on this proc 159 if(usingpetsc.eq.1) then 160 allocate (BDiag(1,1,1)) 161 else 162 allocate (BDiag(nshg,nflow,nflow)) 163 allocate (lhsK(nflow*nflow,nnz_tot)) 164 endif 165 endif 166 if (mod(impl(1),100)/10 .eq. 3) then 167c 168c generate the ebe data fill vectors 169c 170 nedof=nflow*nshape 171 allocate (EGmass(numel,nedof*nedof)) 172 endif 173 174c.......................................... 175 rerr = zero 176 ybar(:,1:ndof) = y(:,1:ndof) 177 do idof=1,5 178 ybar(:,ndof+idof) = y(:,idof)*y(:,idof) 179 enddo 180 ybar(:,ndof+6) = y(:,1)*y(:,2) 181 ybar(:,ndof+7) = y(:,1)*y(:,3) 182 ybar(:,ndof+8) = y(:,2)*y(:,3) 183c......................................... 184 185! change the freestream and inflow eddy viscosity to reflect different 186! levels of freestream turbulence 187! 188! First override boundary condition values 189! USES imapInlet from Mach Control so if that gets conditionaled away 190! it has to know if it is needed here 191! 192 if(isetEV_IC_BC.eq.1) then 193 allocate(vortG(nshg, 4)) 194 call vortGLB(yold, x, shp, shgl, ilwork, vortG) 195 vortmax=maxval(abs(vortG(:,4))) !column 4 is the magnitude of the shear tensor - should actually probably be calld shearmax instead of vortmax 196 write(*,*) "vortmax = ", vortmax 197 198 !Find the maximum shear in the simulation 199 if(numpe.gt.1) then 200 call MPI_ALLREDUCE(vortmax, vortmaxg, 1, 201 & MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ierr ) 202 vortmax = vortmaxg 203 endif 204 205 !Apply eddy viscosity at the inlet 206 do i=1,icountInlet ! for now only coded to catch primary inlet, not blower 207 BC(imapInlet(i),7)=evis_IC_BC 208 enddo 209 210 !Apply eddy viscosity through the quasi-inviscid portion of the domain 211 do i = 1,nshg 212 if(abs(vortG(i,4)).ge.vortmax*0.01) yold(i,6)=evis_IC_BC 213 enddo 214 isclr=1 ! fix scalar 215 call itrBCsclr(yold,ac,iBC,BC,iper,ilwork) 216 217 deallocate(vortG) 218 endif 219c 220c.... loop through the time sequences 221c 222 do 3000 itsq = 1, ntseq 223 itseq = itsq 224c 225c.... set up the current parameters 226c 227 nstp = nstep(itseq) 228 nitr = niter(itseq) 229 LCtime = loctim(itseq) 230c 231 call itrSetup ( y, acold) 232 isclr=0 233 234 niter(itseq)=0 ! count number of flow solves in a step 235 ! (# of iterations) 236 do i=1,seqsize 237 if(stepseq(i).eq.0) niter(itseq)=niter(itseq)+1 238 enddo 239 nitr = niter(itseq) 240c 241c.... determine how many scalar equations we are going to need to solve 242c 243 nsclrsol=nsclr ! total number of scalars solved. At 244 ! some point we probably want to create 245 ! a map, considering stepseq(), to find 246 ! what is actually solved and only 247 ! dimension EGmasst to the appropriate 248 ! size. 249 250 if(nsclrsol.gt.0)allocate (EGmasst(numel*nshape*nshape 251 & ,nsclrsol)) 252 253c 254c.... loop through the time steps 255c 256 ttim(1) = REAL(secs(0.0)) / 100. 257 ttim(2) = secs(0.0) 258 259c tcorecp1 = REAL(secs(0.0)) / 100. 260c tcorewc1 = secs(0.0) 261 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 262 if(myrank.eq.master) then 263 tcorecp1 = TMRC() 264 endif 265 266 rmub=datmat(1,2,1) 267 if(rmutarget.gt.0) then 268 rmue=rmutarget 269 xmulfact=(rmue/rmub)**(1.0/nstp) 270 if(myrank.eq.master) then 271 write(*,*) 'viscosity will by multiplied by ', xmulfact 272 write(*,*) 'to bring it from ', rmub,' down to ', rmue 273 endif 274 datmat(1,2,1)=datmat(1,2,1)/xmulfact ! make first step right 275 else 276 rmue=datmat(1,2,1) ! keep constant 277 xmulfact=one 278 endif 279 if(iramp.eq.1) then 280 call initBCprofileScale(vbc_prof,BC,yold,x) 281! fix the yold values to the reset BC 282 call itrBC (yold, ac, iBC, BC, iper, ilwork) 283 isclr=1 ! fix scalar 284 call itrBCsclr(yold,ac,iBC,BC,iper,ilwork) 285 endif 286 287867 continue 288 289 290c============ Start the loop of time steps============================c 291! edamp2=.99 292! edamp3=0.05 293 deltaInlInv=one/(0.125*0.0254) 294 do 2000 istp = 1, nstp 295 296 297 if(iramp.eq.1) 298 & call BCprofileScale(vbc_prof,BC,yold) 299 300c call rerun_check(stopjob) 301cc if(stopjob.ne.0) goto 2001 302c 303c Decay of scalars 304c 305 if(nsclr.gt.0 .and. tdecay.ne.1) then 306 yold(:,6:ndof)=y(:,6:ndof)*tdecay 307 BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*tdecay 308 endif 309 310 if(nosource.eq.1) BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*0.8 311 312 313c xi=istp*one/nstp 314c datmat(1,2,1)=rmub*(1.0-xi)+xi*rmue 315 datmat(1,2,1)=xmulfact*datmat(1,2,1) 316 317c.... if we have time varying boundary conditions update the values of BC. 318c these will be for time step n+1 so use lstep+1 319c 320 if(itvn.gt.0) call BCint((lstep+1)*Delt(1), shp, shgl, 321 & shpb, shglb, x, BC, iBC) 322 323 if(iLES.gt.0) then 324c 325c.... get dynamic model coefficient 326c 327 ilesmod=iLES/10 328c 329c digit bit set filter rule, 10 bit set model 330c 331 if (ilesmod.eq.0) then ! 0 < iLES < 10 => dyn. model calculated 332 ! at nodes based on discrete filtering 333 call getdmc (yold, shgl, shp, 334 & iper, ilwork, nsons, 335 & ifath, x) 336 endif 337 if (ilesmod .eq. 1) then ! 10 < iLES < 20 => dynamic-mixed 338 ! at nodes based on discrete filtering 339 call bardmc (yold, shgl, shp, 340 & iper, ilwork, 341 & nsons, ifath, x) 342 endif 343 if (ilesmod .eq. 2) then ! 20 < iLES < 30 => dynamic at quad 344 ! pts based on lumped projection filt. 345 call projdmc (yold, shgl, shp, 346 & iper, ilwork, x) 347 endif 348c 349 endif ! endif of iLES 350 351 352c 353c.... set traction BCs for modeled walls 354c 355 if (itwmod.ne.0) then !wallfn check 356 call asbwmod(yold, acold, x, BC, iBC, 357 & iper, ilwork, ifath, velbar) 358 endif 359c 360c.... -----------------------> predictor phase <----------------------- 361c 362 call itrPredict( yold, acold, y, ac ) 363 call itrBC (y, ac, iBC, BC, iper, ilwork) 364 isclr = zero 365 if (nsclr.gt.zero) then 366 do isclr=1,nsclr 367 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 368 enddo 369 endif 370c 371c.... --------------------> multi-corrector phase <-------------------- 372c 373 iter=0 374 ilss=0 ! this is a switch thrown on first solve of LS redistance 375c 376cHACK to make it keep solving RANS until tolerance is reached 377c 378 istop=0 379! blocking extra RANS steps for now iMoreRANS=0 380 iMoreRANS=6 381c 382c find the the RANS portion of the sequence 383c 384 do istepc=1,seqsize 385 if(stepseq(istepc).eq.10) iLastRANS=istepc 386 enddo 387 388 iseqStart=1 3899876 continue 390c 391 do istepc=iseqStart,seqsize 392 icode=stepseq(istepc) 393 if(mod(icode,10).eq.0) then ! this is a solve 394 isolve=icode/10 395 if(isolve.eq.0) then ! flow solve (encoded as 0) 396c 397 etol=epstol(1) 398 iter = iter+1 399 ifuncs(1) = ifuncs(1) + 1 400c 401c.... reset the aerodynamic forces 402c 403 Force(1) = zero 404 Force(2) = zero 405 Force(3) = zero 406 HFlux = zero 407c 408c.... form the element data and solve the matrix problem 409c 410c.... explicit solver 411c 412 if (impl(itseq) .eq. 0) then 413 if (myrank .eq. master) 414 & call error('itrdrv ','impl ',impl(itseq)) 415 endif 416 if (mod(impl(1),100)/10 .eq. 1) then ! sparse solve 417c 418c.... preconditioned sparse matrix GMRES solver 419c 420 lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 421 iprec=lhs 422 nedof = nflow*nshape 423c write(*,*) 'lhs=',lhs 424 if(usingpetsc.eq.1) then 425#if (HAVE_PETSC) 426 call SolGMRp (y, ac, yold, 427 & x, 428 & iBC, BC, 429 & colm, rowp, lhsk, 430 & res, 431 & BDiag, 432 & iper, ilwork, 433 & shp, shgl, 434 & shpb, shglb, solinc, 435 & rerr, fncorp ) 436#else 437 write(*,*) 'exiting because run time input asked for PETSc, not linked in exec' 438 call error('itrdrv ','noPETSc',usingpetsc) 439#endif 440 else 441 call SolGMRs (y, ac, yold, 442 & acold, x, 443 & iBC, BC, 444 & colm, rowp, lhsk, 445 & res, 446 & BDiag, a(mHBrg), a(meBrg), 447 & a(myBrg), a(mRcos), a(mRsin), 448 & iper, ilwork, 449 & shp, shgl, 450 & shpb, shglb, solinc, 451 & rerr) 452 endif 453 else if (mod(impl(1),100)/10 .eq. 2) then ! mfg solve 454c 455c.... preconditioned matrix-free GMRES solver 456c 457 lhs=0 458 iprec = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 459 nedof = 0 460 call SolMFG (y, ac, yold, 461 & acold, x, 462 & iBC, BC, 463 & res, 464 & BDiag, a(mHBrg), a(meBrg), 465 & a(myBrg), a(mRcos), a(mRsin), 466 & iper, ilwork, 467 & shp, shgl, 468 & shpb, shglb, solinc, 469 & rerr) 470c 471 else if (mod(impl(1),100)/10 .eq. 3) then ! ebe solve 472c.... preconditioned ebe matrix GMRES solver 473c 474 475 lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 476 iprec = lhs 477 nedof = nflow*nshape 478c write(*,*) 'lhs=',lhs 479 call SolGMRe (y, ac, yold, 480 & acold, x, 481 & iBC, BC, 482 & EGmass, res, 483 & BDiag, a(mHBrg), a(meBrg), 484 & a(myBrg), a(mRcos), a(mRsin), 485 & iper, ilwork, 486 & shp, shgl, 487 & shpb, shglb, solinc, 488 & rerr) 489 endif 490c 491 else ! solve a scalar (encoded at isclr*10) 492 isclr=isolve 493 etol=epstol(isclr+2) ! note that for both epstol and LHSupd 1 is flow 2 temp isclr+2 for scalars 494 ifuncs(isclr+2) = ifuncs(isclr+2) + 1 495 if((iLSet.eq.2).and.(ilss.eq.0) 496 & .and.(isclr.eq.2)) then 497 ilss=1 ! throw switch (once per step) 498c 499c... copy the first scalar at t=n+1 into the second scalar of the 500c... level sets 501c 502 y(:,6) = yold(:,6) + (y(:,6)-yold(:,6))/alfi 503 y(:,7) = y(:,6) 504 yold(:,7) = y(:,7) 505 ac(:,7) = zero 506c 507 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 508c 509c....store the flow alpha, gamma parameter values and assigm them the 510c....Backward Euler parameters to solve the second levelset scalar 511c 512 alfit=alfi 513 gamit=gami 514 almit=almi 515 alfi = 1 516 gami = 1 517 almi = 1 518 endif 519c 520 lhs = 1 - min(1,mod(ifuncs(isclr+2)-1, 521 & LHSupd(isclr+2))) 522 iprec = lhs 523 istop=0 524 if(usingPETSc.eq.1) then 525#if (HAVE_PETSC) 526 call SolGMRpSclr(y, ac, 527 & x, elDw, 528 & iBC, BC, 529 & colm, rowp, 530 & iper, ilwork, 531 & shp, shgl, 532 & shpb, shglb, rest, 533 & solinc(1,isclr+5),fncorp) 534#else 535 write(*,*) 'exiting because run time input asked for PETSc, not linked in exec' 536 call error('itrdrv ','noPETSc',usingpetsc) 537#endif 538 else 539 call SolGMRSclr(y, ac, yold, 540 & acold, EGmasst(1,isclr), 541 & x, elDw, 542 & iBC, BC, 543 & rest, 544 & a(mHBrg), a(meBrg), 545 & a(myBrg), a(mRcos), a(mRsin), 546 & iper, ilwork, 547 & shp, shgl, 548 & shpb, shglb, solinc(1,isclr+5)) 549 endif 550c 551 endif ! end of scalar type solve 552c 553c 554c.... end of the multi-corrector loop 555c 556 1000 continue !check this 557 558 else ! this is an update (mod did not equal zero) 559 iupdate=icode/10 ! what to update 560 if(iupdate.eq.0) then !update flow 561 call itrCorrect ( y, ac, yold, acold, solinc) 562 call itrBC (y, ac, iBC, BC, iper, ilwork) 563 call tnanq(y, 5, 'y_updbc') 564c Elaine-SPEBC 565 if((irscale.ge.0).and.(myrank.eq.master)) then 566 call genscale(y, x, iBC) 567c call itrBC (y, ac, iBC, BC, iper, ilwork) 568 endif 569 else ! update scalar 570 isclr=iupdate !unless 571 if(iupdate.eq.nsclr+1) isclr=0 572 call itrCorrectSclr ( y, ac, yold, acold, 573 & solinc(1,isclr+5)) 574 if (ilset.eq.2 .and. isclr.eq.2) then 575 fct2=one/almi 576 fct3=one/alfi 577 acold(:,7) = acold(:,7) 578 & + (ac(:,7)-acold(:,7))*fct2 579 yold(:,7) = yold(:,7) 580 & + (y(:,7)-yold(:,7))*fct3 581 call itrBCSclr ( yold, acold, iBC, BC, 582 & iper, ilwork) 583 ac(:,7) = acold(:,7)*(one-almi/gami) 584 y(:,7) = yold(:,7) 585 ac(:,7) = zero 586 if (ivconstraint .eq. 1) then 587 & 588c ... applying the volume constraint 589c 590 call solvecon (y, x, iBC, BC, 591 & iper, ilwork, shp, shgl) 592c 593 endif ! end of volume constraint calculations 594 endif 595 call itrBCSclr ( y, ac, iBC, BC, iper, ilwork) 596 endif 597 endif !end of switch between solve or update 598 enddo ! loop over sequence in step 599 if((istop.lt.0).and.(iMoreRANS.lt.5)) then 600 iMoreRANS=iMoreRANS+1 601 if(myrank.eq.master) write(*,*) 'istop =', istop 602 iseqStart=iLastRANS 603 goto 9876 604 endif 605c 606c Find the solution at the end of the timestep and move it to old 607c 608c.... First to reassign the parameters for the original time integrator scheme 609c 610 if((iLSet.eq.2).and.(ilss.eq.1)) then 611 alfi =alfit 612 gami =gamit 613 almi =almit 614 endif 615 call itrUpdate( yold, acold, y, ac) 616 call itrBC (yold, acold, iBC, BC, iper,ilwork) 617c Elaine-SPEBC 618 if((irscale.ge.0).and.(myrank.eq.master)) then 619 call genscale(yold, x, iBC) 620c call itrBC (y, ac, iBC, BC, iper, ilwork) 621 endif 622 do isclr=1,nsclr 623 call itrBCSclr (yold, acold, iBC, BC, iper, ilwork) 624 enddo 625c 626 istep = istep + 1 627 lstep = lstep + 1 628 ntoutv=max(ntout,100) 629 !boundary flux output moved after the error calculation so 630 !everything can be written out in a single chunk of code - 631 !Nicholas Mati 632 633 !dump TIME SERIES 634 if (exts) then 635 !Write the probe data to disc. Note that IO to disc only 636 !occurs when mod(lstep, nbuff) == 0. However, this 637 !function also does data buffering and must be called 638 !every time step. 639 call TD_bufferData() 640 call TD_writeData(fvarts, .false.) 641 endif 642 643 !Update the Mach Control 644 if(exts .and. exMC) then 645 !Note: the function MC_updateState must be called after 646 !the function TD_bufferData due to dependencies on 647 !vartsbuff(:,:,:). 648 call MC_updateState() 649 call MC_applyBC(BC) 650 call MC_printState() 651 652 !Write the state if a restart is also being written. 653 if(mod(lstep,ntout).eq.0 ) call MC_writeState(lstep) 654 endif 655 656 !update blower control 657 if(BC_enable) then 658 !Update the blower boundary conditions for the next 659 !iteration. 660 call BC_iter(BC) 661 662 !Also write the current phases of the blowers if a 663 !restart is also being written. 664 if(mod(lstep, ntout) == 0) call BC_writePhase(lstep) 665 endif 666 667 !.... Yi Chen Duct geometry8 668 if(isetBlowing_Duct.gt.0)then 669 if(ifixBlowingVel_Duct.eq.0)then 670 if(nstp.gt.nBlowingStepsDuct)then 671 nBlowingStepsDuct = nstp-2 672 endif 673 call setBlowing_Duct2(x,BC,yold,iTurbWall,istp) 674 endif 675 endif 676 !... Yi Chen Duct geometry8 677 678c 679c.... -------------------> error calculation <----------------- 680 if(ierrcalc.eq.1.or.ioybar.eq.1) then 681 tfact=one/istep 682 do idof=1,ndof 683 ybar(:,idof) =tfact*yold(:,idof) + 684 & (one-tfact)*ybar(:,idof) 685 enddo 686c....compute average 687c... ybar(:,ndof+1:ndof+8) is for avg. of square as uu, vv, ww, pp, TT, uv, uw, and vw 688 do idof=1,5 ! avg. of square for 5 flow variables 689 ybar(:,ndof+idof) = tfact*yold(:,idof)**2 + 690 & (one-tfact)*ybar(:,ndof+idof) 691 enddo 692 ybar(:,ndof+6) = tfact*yold(:,1)*yold(:,2) + !uv 693 & (one-tfact)*ybar(:,ndof+6) 694 ybar(:,ndof+7) = tfact*yold(:,1)*yold(:,3) + !uw 695 & (one-tfact)*ybar(:,ndof+7) 696 ybar(:,ndof+8) = tfact*yold(:,2)*yold(:,3) + !vw 697 & (one-tfact)*ybar(:,ndof+8) 698c... compute err 699 rerr(:, 7)=rerr(:, 7)+(yold(:,1)-ybar(:,1))**2 700 rerr(:, 8)=rerr(:, 8)+(yold(:,2)-ybar(:,2))**2 701 rerr(:, 9)=rerr(:, 9)+(yold(:,3)-ybar(:,3))**2 702 rerr(:,10)=rerr(:,10)+(yold(:,4)-ybar(:,4))**2 703 endif 704 705c.. writing ybar field if requested in each restart file 706 707 !here is where we save our averaged field. In some cases we want to 708 !write it less frequently 709 if( (irs >= 1) .and. ( 710 & mod(lstep, ntout) == 0 .or. !Checkpoint 711 & istep == nstp) )then !End of simulation 712 713 if( (mod(lstep, ntoutv) .eq. 0) .and. 714 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 715 & ((nsonmax.eq.1).and.(iLES.gt.0)))) 716 & call rwvelb ('out ', velbar ,ifail) 717 718 !BUG: need to update new_interface to work with SyncIO. 719 !Bflux is presently completely crippled. Note that restar 720 !has also been moved here for readability. 721! call Bflux (yold, acold, x, compute boundary fluxes and print out 722! & shp, shgl, shpb, 723! & shglb, nodflx, ilwork) 724 725 call timer ('Output ') !set up the timer 726 727 !write the solution and time derivative 728 call restar ('out ', yold, acold) 729 730 !Write the distance to wall field in each restart 731 if((istep==nstp) .and. (irans < 0 )) then !d2wall is allocated 732 call write_field(myrank,'a'//char(0),'dwal'//char(0),4, 733 & d2wall,'d'//char(0), nshg, 1, lstep) 734 endif 735 736 !Write the time average in each restart. 737 if(ioybar.eq.1)then 738 call write_field(myrank,'a'//char(0),'ybar'//char(0),4, 739 & ybar,'d'//char(0),nshg,ndof+8,lstep) 740 endif 741 742 !Write the error feild at the end of each step sequence 743 if(ierrcalc.eq.1 .and. istep == nstp) then 744 !smooth the error indicators 745 746 do i=1,ierrsmooth 747 call errsmooth( rerr, x, iper, ilwork, shp, shgl, iBC ) 748 enddo 749 750! call write_error(myrank, lstep, nshg, 10, rerr ) 751 call write_field( 752 & myrank, 'a'//char(0), 'errors'//char(0), 6, 753 & rerr, 'd'//char(0), nshg, 10, lstep) 754 endif 755 756c the following is a future way to have the number of steps in the header...done for posix but not yet for syncio 757c 758c call write_field2(myrank,'a'//char(0),'ybar'//char(0), 759c & 4,ybar,'d'//char(0),nshg,ndof+8, 760c & lstep,istep) 761 endif 762 763 2000 continue !end of NSTEP loop 764 2001 continue 765 766 ttim(1) = REAL(secs(0.0)) / 100. - ttim(1) 767 ttim(2) = secs(0.0) - ttim(2) 768 769c tcorecp2 = REAL(secs(0.0)) / 100. 770c tcorewc2 = secs(0.0) 771 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 772 if(myrank.eq.master) then 773 tcorecp2 = TMRC() 774 write(6,*) 'T(core) cpu = ',tcorecp2-tcorecp1 775 endif 776 777c call wtime 778 779 3000 continue !end of NTSEQ loop 780c 781c.... ----------------------> Post Processing <---------------------- 782c 783c.... print out the last step 784c 785! if ( (irs .ge. 1) .and. ((mod(lstep, ntout) .ne. 0) .or. 786! & (nstp .eq. 0))) then 787! if( (mod(lstep, ntoutv) .eq. 0) .and. 788! & ((irscale.ge.0).or.(itwmod.gt.0) .or. 789! & ((nsonmax.eq.1).and.(iLES.gt.0)))) 790! & call rwvelb ('out ', velbar ,ifail) 791! 792! call Bflux (yold, acold, x, 793! & shp, shgl, shpb, 794! & shglb, nodflx, ilwork) 795! endif 796 797 798 799c if(ioybar.eq.1) then 800c call write_field(myrank,'a'//char(0),'ybar'//char(0),4, 801c & ybar,'d'//char(0),nshg,ndof+8,lstep) 802c endif 803 804c if(iRANS.lt.0 .and. idistcalc.eq.1) then 805c call write_field(myrank,'a'//char(0),'DESd'//char(0),4, 806c & elDw,'d'//char(0),numel,1,lstep) 807c 808c call write_field(myrank,'a'//char(0),'dwal'//char(0),4, 809c & d2wall,'d'//char(0),nshg,1,lstep) 810c endif 811 812c 813c.... close history and aerodynamic forces files 814c 815 if (myrank .eq. master) then 816 close (ihist) 817 close (iforce) 818 819 if(exMC) then 820 call MC_writeState(lstep) 821 call MC_finalize 822 endif 823 824 if(exts) then 825 call TD_writeData(fvarts, .true.) !force the flush of the buffer. 826 call TD_finalize 827 endif 828 endif 829 830 if(BC_enable) then !blower is allocated on all processes. 831 if(mod(lstep, ntout) /= 0) then !May have already written file. 832 call BC_writePhase(lstep) 833 endif 834 835 call BC_finalize 836 endif 837 838 close (iecho) 839 if(iabc==1) deallocate(acs) 840c 841c.... end 842c 843 return 844 end 845 846 847