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 timedata !allows collection of time series 34 use turbSA 35 36 include "common.h" 37 include "mpif.h" 38 include "auxmpi.h" 39 40c 41 dimension y(nshg,ndof), ac(nshg,ndof), 42 & yold(nshg,ndof), acold(nshg,ndof), 43 & x(numnp,nsd), iBC(nshg), 44 & BC(nshg,ndofBC), ilwork(nlwork), 45 & iper(nshg), uold(nshg,nsd) 46c 47 dimension res(nshg,nflow), BDiag(nshg,nflow,nflow), 48 & rest(nshg), solinc(nshg,ndof) 49c 50 dimension shp(MAXTOP,maxsh,MAXQPT), 51 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 52 & shpb(MAXTOP,maxsh,MAXQPT), 53 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 54 real*8 almit, alfit, gamit 55 dimension ifath(numnp), velbar(nfath,ndof), nsons(nfath) 56 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 57 integer, allocatable, dimension(:) :: ivarts 58 integer, allocatable, dimension(:) :: ivartsg 59 real*8, allocatable, dimension(:) :: vartssoln 60 real*8, allocatable, dimension(:) :: vartssolng 61 real*8, allocatable, dimension(:,:,:) :: vartsbuff 62! assuming three profiles to control (inlet, bottom FC and top FC) 63! store velocity profile set via BC 64 real*8 vbc_prof(nshg,3) 65 character*20 fname1, fmt1, fname2, fmt2 66 character*4 fname4c ! 4 characters 67 character*60 fvarts 68 character*5 cname 69 character*10 cname2 70 integer ifuncs(6), iarray(10) 71 72 real*8 elDw(numel) ! element average of DES d variable 73c 74c Here are the data structures for sparse matrix GMRES 75c 76 integer, allocatable, dimension(:,:) :: rowp 77 integer, allocatable, dimension(:) :: colm 78 real*8, allocatable, dimension(:,:) :: lhsK 79 real*8, allocatable, dimension(:,:) :: EGmass 80 real*8, allocatable, dimension(:,:) :: EGmasst 81 82 inquire(file='xyzts.dat',exist=exts) 83 lskeep=lstep 84 if(exts) then 85 86 open(unit=626,file='xyzts.dat',status='old') 87 read(626,*) ntspts, freq, tolpt, iterat, varcod 88 call sTD ! sets data structures 89 do jj=1,ntspts ! read coordinate data where solution desired 90 read(626,*) ptts(jj,1),ptts(jj,2),ptts(jj,3) 91 enddo 92 close(626) 93 94 statptts(:,:) = 0 95 parptts(:,:) = zero 96 varts(:,:) = zero 97 98 allocate (ivarts(ntspts*ndof)) 99 allocate (ivartsg(ntspts*ndof)) 100 allocate (vartssoln(ntspts*ndof)) 101 allocate (vartssolng(ntspts*ndof)) 102 103 nbuff=ntout 104 allocate (vartsbuff(ntspts,ndof,nbuff)) 105 106 if (myrank .eq. master) then 107 do jj=1,ntspts 108 fvarts='varts/varts' 109 fvarts=trim(fvarts)//trim(cname2(jj)) 110 fvarts=trim(fvarts)//trim(cname2(lstep)) 111 fvarts=trim(fvarts)//'.dat' 112 fvarts=trim(fvarts) 113 open(unit=1000+jj, file=fvarts, status='unknown') 114 enddo 115 endif 116 117 endif 118 119c 120c 121c.... open history and aerodynamic forces files 122c 123 if (myrank .eq. master) then 124 open (unit=ihist, file=fhist, status='unknown') 125 open (unit=iforce, file=fforce, status='unknown') 126 endif 127c 128c 129c.... initialize 130c 131 ifuncs = 0 ! func. evaluation counter 132 istep = 0 133 ntotGM = 0 ! number of GMRES iterations 134 time = 0 135 yold = y 136 acold = ac 137 if (mod(impl(1),100)/10 .eq. 1) then 138c 139c generate the sparse data fill vectors 140c 141 allocate (rowp(nshg,nnz)) 142 allocate (colm(nshg+1)) 143 call genadj(colm, rowp, icnt ) ! preprocess the adjacency list 144 145 nnz_tot=icnt ! this is exactly the number of non-zero 146 ! blocks on this proc 147 allocate (lhsK(nflow*nflow,nnz_tot)) 148 endif 149 if (mod(impl(1),100)/10 .eq. 3) then 150c 151c generate the ebe data fill vectors 152c 153 nedof=nflow*nshape 154 allocate (EGmass(numel,nedof*nedof)) 155 endif 156cc 157 158c.......................................... 159 rerr = zero 160 ybar(:,1:ndof) = y(:,1:ndof) 161 do idof=1,5 162 ybar(:,ndof+idof) = y(:,idof)*y(:,idof) 163 enddo 164 ybar(:,ndof+6) = y(:,1)*y(:,2) 165 ybar(:,ndof+7) = y(:,1)*y(:,3) 166 ybar(:,ndof+8) = y(:,2)*y(:,3) 167c......................................... 168 169 170 vbc_prof(:,1:3) = BC(:,3:5) 171 172c 173c.... loop through the time sequences 174c 175 do 3000 itsq = 1, ntseq 176 itseq = itsq 177c 178c.... set up the current parameters 179c 180 nstp = nstep(itseq) 181 nitr = niter(itseq) 182 LCtime = loctim(itseq) 183c 184 call itrSetup ( y, acold) 185 isclr=0 186 187 niter(itseq)=0 ! count number of flow solves in a step 188 ! (# of iterations) 189 do i=1,seqsize 190 if(stepseq(i).eq.0) niter(itseq)=niter(itseq)+1 191 enddo 192 nitr = niter(itseq) 193c 194c.... determine how many scalar equations we are going to need to solve 195c 196 nsclrsol=nsclr ! total number of scalars solved. At 197 ! some point we probably want to create 198 ! a map, considering stepseq(), to find 199 ! what is actually solved and only 200 ! dimension EGmasst to the appropriate 201 ! size. 202 203 if(nsclrsol.gt.0)allocate (EGmasst(numel*nshape*nshape 204 & ,nsclrsol)) 205 206c 207c.... loop through the time steps 208c 209 ttim(1) = REAL(secs(0.0)) / 100. 210 ttim(2) = secs(0.0) 211 212c tcorecp1 = REAL(secs(0.0)) / 100. 213c tcorewc1 = secs(0.0) 214 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 215 if(myrank.eq.0) then 216 tcorecp1 = TMRC() 217 endif 218 219 rmub=datmat(1,2,1) 220 if(rmutarget.gt.0) then 221 rmue=rmutarget 222 xmulfact=(rmue/rmub)**(1.0/nstp) 223 if(myrank.eq.0) then 224 write(*,*) 'viscosity will by multiplied by ', xmulfact 225 write(*,*) 'to bring it from ', rmub,' down to ', rmue 226 endif 227 datmat(1,2,1)=datmat(1,2,1)/xmulfact ! make first step right 228 else 229 rmue=datmat(1,2,1) ! keep constant 230 xmulfact=one 231 endif 232 if(iramp.eq.1) then 233 call initBCprofileScale(vbc_prof,BC,yold,x) 234! fix the yold values to the reset BC 235 call itrBC (yold, ac, iBC, BC, iper, ilwork) 236 isclr=1 ! fix scalar 237 call itrBCsclr(yold,ac,iBC,BC,iper,ilwork) 238 endif 239 240867 continue 241 do 2000 istp = 1, nstp 242 if(iramp.eq.1) 243 & call BCprofileScale(vbc_prof,BC,yold) 244 245c call rerun_check(stopjob) 246cc if(stopjob.ne.0) goto 2001 247c 248c Decay of scalars 249c 250 if(nsclr.gt.0 .and. tdecay.ne.1) then 251 yold(:,6:ndof)=y(:,6:ndof)*tdecay 252 BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*tdecay 253 endif 254 255 if(nosource.eq.1) BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*0.8 256 257 258c xi=istp*one/nstp 259c datmat(1,2,1)=rmub*(1.0-xi)+xi*rmue 260 datmat(1,2,1)=xmulfact*datmat(1,2,1) 261 262c.... if we have time varying boundary conditions update the values of BC. 263c these will be for time step n+1 so use lstep+1 264c 265 if(itvn.gt.0) call BCint((lstep+1)*Delt(1), shp, shgl, 266 & shpb, shglb, x, BC, iBC) 267 268 if(iLES.gt.0) then 269c 270c.... get dynamic model coefficient 271c 272 ilesmod=iLES/10 273c 274c digit bit set filter rule, 10 bit set model 275c 276 if (ilesmod.eq.0) then ! 0 < iLES < 10 => dyn. model calculated 277 ! at nodes based on discrete filtering 278 call getdmc (yold, shgl, shp, 279 & iper, ilwork, nsons, 280 & ifath, x) 281 endif 282 if (ilesmod .eq. 1) then ! 10 < iLES < 20 => dynamic-mixed 283 ! at nodes based on discrete filtering 284 call bardmc (yold, shgl, shp, 285 & iper, ilwork, 286 & nsons, ifath, x) 287 endif 288 if (ilesmod .eq. 2) then ! 20 < iLES < 30 => dynamic at quad 289 ! pts based on lumped projection filt. 290 call projdmc (yold, shgl, shp, 291 & iper, ilwork, x) 292 endif 293c 294 endif 295 296 297c 298c.... set traction BCs for modeled walls 299c 300 if (itwmod.ne.0) then !wallfn check 301 call asbwmod(yold, acold, x, BC, iBC, 302 & iper, ilwork, ifath, velbar) 303 endif 304c 305c.... -----------------------> predictor phase <----------------------- 306c 307 call itrPredict( yold, acold, y, ac ) 308 call itrBC (y, ac, iBC, BC, iper, ilwork) 309 isclr = zero 310 if (nsclr.gt.zero) then 311 do isclr=1,nsclr 312 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 313 enddo 314 endif 315c 316c.... --------------------> multi-corrector phase <-------------------- 317c 318 iter=0 319 ilss=0 ! this is a switch thrown on first solve of LS redistance 320c 321cHACK to make it keep solving RANS until tolerance is reached 322c 323 istop=0 324 iMoreRANS=0 325c 326c find the the RANS portion of the sequence 327c 328 do istepc=1,seqsize 329 if(stepseq(istepc).eq.10) iLastRANS=istepc 330 enddo 331 332 iseqStart= 1 3339876 continue 334c 335 do istepc=iseqStart,seqsize 336 icode=stepseq(istepc) 337 if(mod(icode,10).eq.0) then ! this is a solve 338 isolve=icode/10 339 if(isolve.eq.0) then ! flow solve (encoded as 0) 340c 341 etol=epstol(1) 342 iter = iter+1 343 ifuncs(1) = ifuncs(1) + 1 344c 345c.... reset the aerodynamic forces 346c 347 Force(1) = zero 348 Force(2) = zero 349 Force(3) = zero 350 HFlux = zero 351c 352c.... form the element data and solve the matrix problem 353c 354c.... explicit solver 355c 356 if (impl(itseq) .eq. 0) then 357 if (myrank .eq. master) 358 & call error('itrdrv ','impl ',impl(itseq)) 359 endif 360 if (mod(impl(1),100)/10 .eq. 1) then ! sparse solve 361c 362c.... preconditioned sparse matrix GMRES solver 363c 364 lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 365 iprec=lhs 366 nedof = nflow*nshape 367c write(*,*) 'lhs=',lhs 368 call SolGMRs (y, ac, yold, 369 & acold, x, 370 & iBC, BC, 371 & colm, rowp, lhsk, 372 & res, 373 & BDiag, a(mHBrg), a(meBrg), 374 & a(myBrg), a(mRcos), a(mRsin), 375 & iper, ilwork, 376 & shp, shgl, 377 & shpb, shglb, solinc, 378 & rerr) 379 else if (mod(impl(1),100)/10 .eq. 2) then ! mfg solve 380c 381c.... preconditioned matrix-free GMRES solver 382c 383 lhs=0 384 iprec = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 385 nedof = 0 386 call SolMFG (y, ac, yold, 387 & acold, x, 388 & iBC, BC, 389 & res, 390 & BDiag, a(mHBrg), a(meBrg), 391 & a(myBrg), a(mRcos), a(mRsin), 392 & iper, ilwork, 393 & shp, shgl, 394 & shpb, shglb, solinc, 395 & rerr) 396c 397 else if (mod(impl(1),100)/10 .eq. 3) then ! ebe solve 398c.... preconditioned ebe matrix GMRES solver 399c 400 401 lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 402 iprec = lhs 403 nedof = nflow*nshape 404c write(*,*) 'lhs=',lhs 405 call SolGMRe (y, ac, yold, 406 & acold, x, 407 & iBC, BC, 408 & EGmass, res, 409 & BDiag, a(mHBrg), a(meBrg), 410 & a(myBrg), a(mRcos), a(mRsin), 411 & iper, ilwork, 412 & shp, shgl, 413 & shpb, shglb, solinc, 414 & rerr) 415 endif 416c 417 else ! solve a scalar (encoded at isclr*10) 418 ifuncs(isclr+2) = ifuncs(isclr+2) + 1 419 etol=epstol(isclr+1) 420 isclr=isolve 421 if((iLSet.eq.2).and.(ilss.eq.0) 422 & .and.(isclr.eq.2)) then 423 ilss=1 ! throw switch (once per step) 424c 425c... copy the first scalar at t=n+1 into the second scalar of the 426c... level sets 427c 428 y(:,6) = yold(:,6) + (y(:,6)-yold(:,6))/alfi 429 y(:,7) = y(:,6) 430 yold(:,7) = y(:,7) 431 ac(:,7) = zero 432c 433 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 434c 435c....store the flow alpha, gamma parameter values and assigm them the 436c....Backward Euler parameters to solve the second levelset scalar 437c 438 alfit=alfi 439 gamit=gami 440 almit=almi 441 alfi = 1 442 gami = 1 443 almi = 1 444 endif 445c 446 lhs = 1 - min(1,mod(ifuncs(isclr+2)-1, 447 & LHSupd(isclr+2))) 448 iprec = lhs 449 istop=0 450 call SolGMRSclr(y, ac, yold, 451 & acold, EGmasst(1,isclr), 452 & x, elDw, 453 & iBC, BC, 454 & rest, 455 & a(mHBrg), a(meBrg), 456 & a(myBrg), a(mRcos), a(mRsin), 457 & iper, ilwork, 458 & shp, shgl, 459 & shpb, shglb, solinc(1,isclr+5)) 460c 461 endif ! end of scalar type solve 462c 463c 464c.... end of the multi-corrector loop 465c 466 1000 continue !check this 467 468 else ! this is an update (mod did not equal zero) 469 iupdate=icode/10 ! what to update 470 if(iupdate.eq.0) then !update flow 471 call itrCorrect ( y, ac, yold, acold, solinc) 472 call itrBC (y, ac, iBC, BC, iper, ilwork) 473 call tnanq(y, 5, 'y_updbc') 474c Elaine-SPEBC 475 if((irscale.ge.0).and.(myrank.eq.master)) then 476 call genscale(y, x, iBC) 477c call itrBC (y, ac, iBC, BC, iper, ilwork) 478 endif 479 else ! update scalar 480 isclr=iupdate !unless 481 if(iupdate.eq.nsclr+1) isclr=0 482 call itrCorrectSclr ( y, ac, yold, acold, 483 & solinc(1,isclr+5)) 484 if (ilset.eq.2 .and. isclr.eq.2) then 485 fct2=one/almi 486 fct3=one/alfi 487 acold(:,7) = acold(:,7) 488 & + (ac(:,7)-acold(:,7))*fct2 489 yold(:,7) = yold(:,7) 490 & + (y(:,7)-yold(:,7))*fct3 491 call itrBCSclr ( yold, acold, iBC, BC, 492 & iper, ilwork) 493 ac(:,7) = acold(:,7)*(one-almi/gami) 494 y(:,7) = yold(:,7) 495 ac(:,7) = zero 496 if (ivconstraint .eq. 1) then 497 & 498c ... applying the volume constraint 499c 500 call solvecon (y, x, iBC, BC, 501 & iper, ilwork, shp, shgl) 502c 503 endif ! end of volume constraint calculations 504 endif 505 call itrBCSclr ( y, ac, iBC, BC, iper, ilwork) 506 endif 507 endif !end of switch between solve or update 508 enddo ! loop over sequence in step 509 if((istop.lt.0).and.(iMoreRANS.lt.5)) then 510 iMoreRANS=iMoreRANS+1 511 if(myrank.eq.0) write(*,*) 'istop =', istop 512 iseqStart=iLastRANS 513 goto 9876 514 endif 515c 516c Find the solution at the end of the timestep and move it to old 517c 518c.... First to reassign the parameters for the original time integrator scheme 519c 520 if((iLSet.eq.2).and.(ilss.eq.1)) then 521 alfi =alfit 522 gami =gamit 523 almi =almit 524 endif 525 call itrUpdate( yold, acold, y, ac) 526 call itrBC (yold, acold, iBC, BC, iper,ilwork) 527c Elaine-SPEBC 528 if((irscale.ge.0).and.(myrank.eq.master)) then 529 call genscale(yold, x, iBC) 530c call itrBC (y, ac, iBC, BC, iper, ilwork) 531 endif 532 do isclr=1,nsclr 533 call itrBCSclr (yold, acold, iBC, BC, iper, ilwork) 534 enddo 535c 536 istep = istep + 1 537c 538c.... compute boundary fluxes and print out 539c 540 lstep = lstep + 1 541 ntoutv=max(ntout,100) 542 if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then 543c 544c here is where we save our averaged field. In some cases we want to 545c write it less frequently 546 if( (mod(lstep, ntoutv) .eq. 0) .and. 547 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 548 & ((nsonmax.eq.1).and.(iLES.gt.0)))) 549 & call rwvelb ('out ', velbar ,ifail) 550 551 call Bflux (yold, acold, x, 552 & shp, shgl, shpb, 553 & shglb, nodflx, ilwork) 554 endif 555c 556c.... end of the NSTEP and NTSEQ loops 557c 558 559c... dump TIME SERIES 560 561 if (exts) then 562 if (mod(lstep-1,freq).eq.0) then 563 564 if (numpe > 1) then 565 do jj = 1, ntspts 566 vartssoln((jj-1)*ndof+1:jj*ndof)=varts(jj,:) 567 ivarts=zero 568 enddo 569 do k=1,ndof*ntspts 570 if(vartssoln(k).ne.zero) ivarts(k)=1 571 enddo 572 call MPI_REDUCE(vartssoln, vartssolng, ndof*ntspts, 573 & MPI_DOUBLE_PRECISION, MPI_SUM, master, 574 & MPI_COMM_WORLD, ierr) 575 576 call MPI_REDUCE(ivarts, ivartsg, ndof*ntspts, 577 & MPI_INTEGER, MPI_SUM, master, 578 & MPI_COMM_WORLD, ierr) 579 580 if (myrank.eq.zero) then 581 do jj = 1, ntspts 582 583 indxvarts = (jj-1)*ndof 584 do k=1,ndof 585 if(ivartsg(indxvarts+k).ne.0) then ! none of the vartssoln(parts) were non zero 586 varts(jj,k)=vartssolng(indxvarts+k)/ 587 & ivartsg(indxvarts+k) 588 endif 589 enddo 590 enddo 591 endif !only on master 592 endif !only if numpe > 1 593 594 if( istp.eq. nstp) then !make sure incomplete buffers get purged. 595 icheck=mod(nstp,nbuff) 596 if(icheck.ne.0) nbuff=icheck 597 endif 598 599 if (myrank.eq.zero) then 600 k=mod(lstep,nbuff) 601 if(k.eq.0) k=nbuff 602 do jj = 1, ntspts 603 vartsbuff(jj,1:5,k)=varts(jj,1:5) 604 enddo 605 if(k.eq. nbuff) then 606 do jj = 1, ntspts 607 ifile = 1000+jj 608 do ibuf=1,nbuff 609 write(ifile,555) lstep-1 -nbuff+ibuf, 610 & (vartsbuff(jj,k,ibuf), k=1,5) ! assuming ndof to be 5 611 enddo ! buff empty 612 613c call flush(ifile) 614 enddo ! jj ntspts 615 endif !only dump when buffer full 616 endif !only on master 617 618 varts(:,:) = zero ! reset the array for next step 619 620 555 format(i6,5(2x,E12.5e2)) 621 622 endif 623 endif 624 625c 626c.... -------------------> error calculation <----------------- 627 if(ierrcalc.eq.1.or.ioybar.eq.1) then 628 tfact=one/istep 629 do idof=1,ndof 630 ybar(:,idof) =tfact*yold(:,idof) + 631 & (one-tfact)*ybar(:,idof) 632 enddo 633c....compute average 634c... ybar(:,ndof+1:ndof+8) is for avg. of square as uu, vv, ww, pp, TT, uv, uw, and vw 635 do idof=1,5 ! avg. of square for 5 flow variables 636 ybar(:,ndof+idof) = tfact*yold(:,idof)**2 + 637 & (one-tfact)*ybar(:,ndof+idof) 638 enddo 639 ybar(:,ndof+6) = tfact*yold(:,1)*yold(:,2) + !uv 640 & (one-tfact)*ybar(:,ndof+6) 641 ybar(:,ndof+7) = tfact*yold(:,1)*yold(:,3) + !uw 642 & (one-tfact)*ybar(:,ndof+7) 643 ybar(:,ndof+8) = tfact*yold(:,2)*yold(:,3) + !vw 644 & (one-tfact)*ybar(:,ndof+8) 645c... compute err 646c rerr(:, 7)=rerr(:, 7)+(yold(:,1)-ybar(:,1))**2 647c rerr(:, 8)=rerr(:, 8)+(yold(:,2)-ybar(:,2))**2 648c rerr(:, 9)=rerr(:, 9)+(yold(:,3)-ybar(:,3))**2 649c rerr(:,10)=rerr(:,10)+(yold(:,4)-ybar(:,4))**2 650 endif 651c........................................................................... 652 2000 continue 653 2001 continue 654 655 ttim(1) = REAL(secs(0.0)) / 100. - ttim(1) 656 ttim(2) = secs(0.0) - ttim(2) 657 658c tcorecp2 = REAL(secs(0.0)) / 100. 659c tcorewc2 = secs(0.0) 660 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 661 if(myrank.eq.0) then 662 tcorecp2 = TMRC() 663 write(6,*) 'T(core) cpu = ',tcorecp2-tcorecp1 664 endif 665 666c call wtime 667 668 3000 continue 669c 670c.... ----------------------> Post Processing <---------------------- 671c 672c.... print out the last step 673c 674 if ((irs .ge. 1) .and. ((mod(lstep, ntout) .ne. 0) .or. 675 & (nstp .eq. 0))) then 676 if( (mod(lstep, ntoutv) .eq. 0) .and. 677 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 678 & ((nsonmax.eq.1).and.(iLES.gt.0)))) 679 & call rwvelb ('out ', velbar ,ifail) 680 681 call Bflux (yold, acold, x, 682 & shp, shgl, shpb, 683 & shglb, nodflx, ilwork) 684 endif 685c 686 if(ierrcalc.eq.1) then 687c 688c.....smooth the error indicators 689c 690c ! errsmooth is currently available only for incompressible code 691c do i=1,ierrsmooth 692c call errsmooth( rerr, x, iper, ilwork, shp, shgl, iBC ) 693c end do 694 695 call write_error(myrank, lstep, nshg, 10, rerr ) 696 endif 697 698 if(ioybar.eq.1) then 699 call write_field(myrank,'a','ybar',4, 700 & ybar,'d',nshg,ndof+8,lstep) 701 endif 702 703 if(iRANS.lt.0 .and. idistcalc.eq.1) then 704 call write_field(myrank,'a','DESd',4, 705 & elDw,'d',numel,1,lstep) 706 707 call write_field(myrank,'a','dwal',4, 708 & d2wall,'d',nshg,1,lstep) 709 endif 710 711c 712c.... close history and aerodynamic forces files 713c 714 if (myrank .eq. master) then 715 close (ihist) 716 close (iforce) 717 if(exts) then 718 deallocate(ivarts) 719 deallocate(ivartsg) 720 deallocate(vartssoln) 721 deallocate(vartssolng) 722 do jj=1,ntspts 723 close(1000+jj) 724 enddo 725 endif 726 endif 727 close (iecho) 728 if(iabc==1) deallocate(acs) 729c 730c.... end 731c 732 return 733 end 734