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