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