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