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