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 timedata !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 call SolGMRp (y, ac, yold, 426 & x, 427 & iBC, BC, 428 & colm, rowp, lhsk, 429 & res, 430 & BDiag, 431 & iper, ilwork, 432 & shp, shgl, 433 & shpb, shglb, solinc, 434 & rerr, fncorp ) 435 else 436 call SolGMRs (y, ac, yold, 437 & acold, x, 438 & iBC, BC, 439 & colm, rowp, lhsk, 440 & res, 441 & BDiag, a(mHBrg), a(meBrg), 442 & a(myBrg), a(mRcos), a(mRsin), 443 & iper, ilwork, 444 & shp, shgl, 445 & shpb, shglb, solinc, 446 & rerr) 447 endif 448 else if (mod(impl(1),100)/10 .eq. 2) then ! mfg solve 449c 450c.... preconditioned matrix-free GMRES solver 451c 452 lhs=0 453 iprec = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 454 nedof = 0 455 call SolMFG (y, ac, yold, 456 & acold, x, 457 & iBC, BC, 458 & res, 459 & BDiag, a(mHBrg), a(meBrg), 460 & a(myBrg), a(mRcos), a(mRsin), 461 & iper, ilwork, 462 & shp, shgl, 463 & shpb, shglb, solinc, 464 & rerr) 465c 466 else if (mod(impl(1),100)/10 .eq. 3) then ! ebe solve 467c.... preconditioned ebe matrix GMRES solver 468c 469 470 lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 471 iprec = lhs 472 nedof = nflow*nshape 473c write(*,*) 'lhs=',lhs 474 call SolGMRe (y, ac, yold, 475 & acold, x, 476 & iBC, BC, 477 & EGmass, res, 478 & BDiag, a(mHBrg), a(meBrg), 479 & a(myBrg), a(mRcos), a(mRsin), 480 & iper, ilwork, 481 & shp, shgl, 482 & shpb, shglb, solinc, 483 & rerr) 484 endif 485c 486 else ! solve a scalar (encoded at isclr*10) 487 isclr=isolve 488 etol=epstol(isclr+2) ! note that for both epstol and LHSupd 1 is flow 2 temp isclr+2 for scalars 489 ifuncs(isclr+2) = ifuncs(isclr+2) + 1 490 if((iLSet.eq.2).and.(ilss.eq.0) 491 & .and.(isclr.eq.2)) then 492 ilss=1 ! throw switch (once per step) 493c 494c... copy the first scalar at t=n+1 into the second scalar of the 495c... level sets 496c 497 y(:,6) = yold(:,6) + (y(:,6)-yold(:,6))/alfi 498 y(:,7) = y(:,6) 499 yold(:,7) = y(:,7) 500 ac(:,7) = zero 501c 502 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 503c 504c....store the flow alpha, gamma parameter values and assigm them the 505c....Backward Euler parameters to solve the second levelset scalar 506c 507 alfit=alfi 508 gamit=gami 509 almit=almi 510 alfi = 1 511 gami = 1 512 almi = 1 513 endif 514c 515 lhs = 1 - min(1,mod(ifuncs(isclr+2)-1, 516 & LHSupd(isclr+2))) 517 iprec = lhs 518 istop=0 519 if(usingPETSc.eq.1) then 520 call SolGMRpSclr(y, ac, 521 & x, elDw, 522 & iBC, BC, 523 & colm, rowp, 524 & iper, ilwork, 525 & shp, shgl, 526 & shpb, shglb, rest, 527 & solinc(1,isclr+5),fncorp) 528 529 else 530 call SolGMRSclr(y, ac, yold, 531 & acold, EGmasst(1,isclr), 532 & x, elDw, 533 & iBC, BC, 534 & rest, 535 & a(mHBrg), a(meBrg), 536 & a(myBrg), a(mRcos), a(mRsin), 537 & iper, ilwork, 538 & shp, shgl, 539 & shpb, shglb, solinc(1,isclr+5)) 540 endif 541c 542 endif ! end of scalar type solve 543c 544c 545c.... end of the multi-corrector loop 546c 547 1000 continue !check this 548 549 else ! this is an update (mod did not equal zero) 550 iupdate=icode/10 ! what to update 551 if(iupdate.eq.0) then !update flow 552 call itrCorrect ( y, ac, yold, acold, solinc) 553 call itrBC (y, ac, iBC, BC, iper, ilwork) 554 call tnanq(y, 5, 'y_updbc') 555c Elaine-SPEBC 556 if((irscale.ge.0).and.(myrank.eq.master)) then 557 call genscale(y, x, iBC) 558c call itrBC (y, ac, iBC, BC, iper, ilwork) 559 endif 560 else ! update scalar 561 isclr=iupdate !unless 562 if(iupdate.eq.nsclr+1) isclr=0 563 call itrCorrectSclr ( y, ac, yold, acold, 564 & solinc(1,isclr+5)) 565 if (ilset.eq.2 .and. isclr.eq.2) then 566 fct2=one/almi 567 fct3=one/alfi 568 acold(:,7) = acold(:,7) 569 & + (ac(:,7)-acold(:,7))*fct2 570 yold(:,7) = yold(:,7) 571 & + (y(:,7)-yold(:,7))*fct3 572 call itrBCSclr ( yold, acold, iBC, BC, 573 & iper, ilwork) 574 ac(:,7) = acold(:,7)*(one-almi/gami) 575 y(:,7) = yold(:,7) 576 ac(:,7) = zero 577 if (ivconstraint .eq. 1) then 578 & 579c ... applying the volume constraint 580c 581 call solvecon (y, x, iBC, BC, 582 & iper, ilwork, shp, shgl) 583c 584 endif ! end of volume constraint calculations 585 endif 586 call itrBCSclr ( y, ac, iBC, BC, iper, ilwork) 587 endif 588 endif !end of switch between solve or update 589 enddo ! loop over sequence in step 590 if((istop.lt.0).and.(iMoreRANS.lt.5)) then 591 iMoreRANS=iMoreRANS+1 592 if(myrank.eq.master) write(*,*) 'istop =', istop 593 iseqStart=iLastRANS 594 goto 9876 595 endif 596c 597c Find the solution at the end of the timestep and move it to old 598c 599c.... First to reassign the parameters for the original time integrator scheme 600c 601 if((iLSet.eq.2).and.(ilss.eq.1)) then 602 alfi =alfit 603 gami =gamit 604 almi =almit 605 endif 606 call itrUpdate( yold, acold, y, ac) 607 call itrBC (yold, acold, iBC, BC, iper,ilwork) 608c Elaine-SPEBC 609 if((irscale.ge.0).and.(myrank.eq.master)) then 610 call genscale(yold, x, iBC) 611c call itrBC (y, ac, iBC, BC, iper, ilwork) 612 endif 613 do isclr=1,nsclr 614 call itrBCSclr (yold, acold, iBC, BC, iper, ilwork) 615 enddo 616c 617 istep = istep + 1 618 lstep = lstep + 1 619 ntoutv=max(ntout,100) 620 !boundary flux output moved after the error calculation so 621 !everything can be written out in a single chunk of code - 622 !Nicholas Mati 623 624 !dump TIME SERIES 625 if (exts) then 626 !Write the probe data to disc. Note that IO to disc only 627 !occurs when mod(lstep, nbuff) == 0. However, this 628 !function also does data buffering and must be called 629 !every time step. 630 call TD_bufferData() 631 call TD_writeData(fvarts, .false.) 632 endif 633 634 !Update the Mach Control 635 if(exts .and. exMC) then 636 !Note: the function MC_updateState must be called after 637 !the function TD_bufferData due to dependencies on 638 !vartsbuff(:,:,:). 639 call MC_updateState() 640 call MC_applyBC(BC) 641 call MC_printState() 642 643 !Write the state if a restart is also being written. 644 if(mod(lstep,ntout).eq.0 ) call MC_writeState(lstep) 645 endif 646 647 !update blower control 648 if(BC_enable) then 649 !Update the blower boundary conditions for the next 650 !iteration. 651 call BC_iter(BC) 652 653 !Also write the current phases of the blowers if a 654 !restart is also being written. 655 if(mod(lstep, ntout) == 0) call BC_writePhase(lstep) 656 endif 657 658 !.... Yi Chen Duct geometry8 659 if(isetBlowing_Duct.gt.0)then 660 if(ifixBlowingVel_Duct.eq.0)then 661 if(nstp.gt.nBlowingStepsDuct)then 662 nBlowingStepsDuct = nstp-2 663 endif 664 call setBlowing_Duct2(x,BC,yold,iTurbWall,istp) 665 endif 666 endif 667 !... Yi Chen Duct geometry8 668 669c 670c.... -------------------> error calculation <----------------- 671 if(ierrcalc.eq.1.or.ioybar.eq.1) then 672 tfact=one/istep 673 do idof=1,ndof 674 ybar(:,idof) =tfact*yold(:,idof) + 675 & (one-tfact)*ybar(:,idof) 676 enddo 677c....compute average 678c... ybar(:,ndof+1:ndof+8) is for avg. of square as uu, vv, ww, pp, TT, uv, uw, and vw 679 do idof=1,5 ! avg. of square for 5 flow variables 680 ybar(:,ndof+idof) = tfact*yold(:,idof)**2 + 681 & (one-tfact)*ybar(:,ndof+idof) 682 enddo 683 ybar(:,ndof+6) = tfact*yold(:,1)*yold(:,2) + !uv 684 & (one-tfact)*ybar(:,ndof+6) 685 ybar(:,ndof+7) = tfact*yold(:,1)*yold(:,3) + !uw 686 & (one-tfact)*ybar(:,ndof+7) 687 ybar(:,ndof+8) = tfact*yold(:,2)*yold(:,3) + !vw 688 & (one-tfact)*ybar(:,ndof+8) 689c... compute err 690 rerr(:, 7)=rerr(:, 7)+(yold(:,1)-ybar(:,1))**2 691 rerr(:, 8)=rerr(:, 8)+(yold(:,2)-ybar(:,2))**2 692 rerr(:, 9)=rerr(:, 9)+(yold(:,3)-ybar(:,3))**2 693 rerr(:,10)=rerr(:,10)+(yold(:,4)-ybar(:,4))**2 694 endif 695 696c.. writing ybar field if requested in each restart file 697 698 !here is where we save our averaged field. In some cases we want to 699 !write it less frequently 700 if( (irs >= 1) .and. ( 701 & mod(lstep, ntout) == 0 .or. !Checkpoint 702 & istep == nstp) )then !End of simulation 703 704 if( (mod(lstep, ntoutv) .eq. 0) .and. 705 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 706 & ((nsonmax.eq.1).and.(iLES.gt.0)))) 707 & call rwvelb ('out ', velbar ,ifail) 708 709 !BUG: need to update new_interface to work with SyncIO. 710 !Bflux is presently completely crippled. Note that restar 711 !has also been moved here for readability. 712! call Bflux (yold, acold, x, compute boundary fluxes and print out 713! & shp, shgl, shpb, 714! & shglb, nodflx, ilwork) 715 716 call timer ('Output ') !set up the timer 717 718 !write the solution and time derivative 719 call restar ('out ', yold, acold) 720 721 !Write the distance to wall field in each restart 722 if((istep==nstp) .and. (irans < 0 )) then !d2wall is allocated 723 call write_field(myrank,'a'//char(0),'dwal'//char(0),4, 724 & d2wall,'d'//char(0), nshg, 1, lstep) 725 endif 726 727 !Write the time average in each restart. 728 if(ioybar.eq.1)then 729 call write_field(myrank,'a'//char(0),'ybar'//char(0),4, 730 & ybar,'d'//char(0),nshg,ndof+8,lstep) 731 endif 732 733 !Write the error feild at the end of each step sequence 734 if(ierrcalc.eq.1 .and. istep == nstp) then 735 !smooth the error indicators 736 737 do i=1,ierrsmooth 738 call errsmooth( rerr, x, iper, ilwork, shp, shgl, iBC ) 739 enddo 740 741! call write_error(myrank, lstep, nshg, 10, rerr ) 742 call write_field( 743 & myrank, 'a'//char(0), 'errors'//char(0), 6, 744 & rerr, 'd'//char(0), nshg, 10, lstep) 745 endif 746 747c the following is a future way to have the number of steps in the header...done for posix but not yet for syncio 748c 749c call write_field2(myrank,'a'//char(0),'ybar'//char(0), 750c & 4,ybar,'d'//char(0),nshg,ndof+8, 751c & lstep,istep) 752 endif 753 754 2000 continue !end of NSTEP loop 755 2001 continue 756 757 ttim(1) = REAL(secs(0.0)) / 100. - ttim(1) 758 ttim(2) = secs(0.0) - ttim(2) 759 760c tcorecp2 = REAL(secs(0.0)) / 100. 761c tcorewc2 = secs(0.0) 762 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 763 if(myrank.eq.master) then 764 tcorecp2 = TMRC() 765 write(6,*) 'T(core) cpu = ',tcorecp2-tcorecp1 766 endif 767 768c call wtime 769 770 3000 continue !end of NTSEQ loop 771c 772c.... ----------------------> Post Processing <---------------------- 773c 774c.... print out the last step 775c 776! if ( (irs .ge. 1) .and. ((mod(lstep, ntout) .ne. 0) .or. 777! & (nstp .eq. 0))) then 778! if( (mod(lstep, ntoutv) .eq. 0) .and. 779! & ((irscale.ge.0).or.(itwmod.gt.0) .or. 780! & ((nsonmax.eq.1).and.(iLES.gt.0)))) 781! & call rwvelb ('out ', velbar ,ifail) 782! 783! call Bflux (yold, acold, x, 784! & shp, shgl, shpb, 785! & shglb, nodflx, ilwork) 786! endif 787 788 789 790c if(ioybar.eq.1) then 791c call write_field(myrank,'a'//char(0),'ybar'//char(0),4, 792c & ybar,'d'//char(0),nshg,ndof+8,lstep) 793c endif 794 795c if(iRANS.lt.0 .and. idistcalc.eq.1) then 796c call write_field(myrank,'a'//char(0),'DESd'//char(0),4, 797c & elDw,'d'//char(0),numel,1,lstep) 798c 799c call write_field(myrank,'a'//char(0),'dwal'//char(0),4, 800c & d2wall,'d'//char(0),nshg,1,lstep) 801c endif 802 803c 804c.... close history and aerodynamic forces files 805c 806 if (myrank .eq. master) then 807 close (ihist) 808 close (iforce) 809 810 if(exMC) then 811 call MC_writeState(lstep) 812 call MC_finalize 813 endif 814 815 if(exts) then 816 call TD_writeData(fvarts, .true.) !force the flush of the buffer. 817 call TD_finalize 818 endif 819 endif 820 821 if(BC_enable) then !blower is allocated on all processes. 822 if(mod(lstep, ntout) /= 0) then !May have already written file. 823 call BC_writePhase(lstep) 824 endif 825 826 call BC_finalize 827 endif 828 829 close (iecho) 830 if(iabc==1) deallocate(acs) 831c 832c.... end 833c 834 return 835 end 836 837 838