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