1 subroutine itrdrv (y, ac, 2 & uold, x, 3 & iBC, BC, 4 & iper, ilwork, shp, 5 & shgl, shpb, shglb, 6 & ifath, velbar, nsons ) 7c 8c---------------------------------------------------------------------- 9c 10c This iterative driver is the semi-discrete, predictor multi-corrector 11c algorithm. It contains the Hulbert Generalized Alpha method which 12c is 2nd order accurate for Rho_inf from 0 to 1. The method can be 13c made first-order accurate by setting Rho_inf=-1. It uses CGP and 14c GMRES iterative solvers. 15c 16c working arrays: 17c y (nshg,ndof) : Y variables 18c x (nshg,nsd) : node coordinates 19c iBC (nshg) : BC codes 20c BC (nshg,ndofBC) : BC constraint parameters 21c iper (nshg) : periodicity table 22c 23c 24c Zdenek Johan, Winter 1991. (Fortran 90) 25c Alberto Figueroa, Winter 2004. CMM-FSI 26c Irene Vignon, Fall 2004. Impedance BC 27c---------------------------------------------------------------------- 28c 29 use pvsQbi !gives us splag (the spmass at the end of this run 30 use specialBC !gives us itvn 31 use timedata !allows collection of time series 32 use convolImpFlow !for Imp bc 33 use convolRCRFlow !for RCR bc 34!MR CHANGE 35 use turbsa ! used to access d2wall 36!MR CHANGE END 37 use iso_c_binding 38 39c use readarrays !reads in uold and acold 40 41 include "common.h" 42 include "mpif.h" 43 include "auxmpi.h" 44c 45 46 47 real*8 y(nshg,ndof), ac(nshg,ndof), 48 & yold(nshg,ndof), acold(nshg,ndof), 49 & u(nshg,nsd), uold(nshg,nsd), 50 & x(numnp,nsd), solinc(nshg,ndof), 51 & BC(nshg,ndofBC), tf(nshg,ndof), 52 & GradV(nshg,nsdsq) 53 54c 55 real*8 res(nshg,ndof) 56c 57 real*8 shp(MAXTOP,maxsh,MAXQPT), 58 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 59 & shpb(MAXTOP,maxsh,MAXQPT), 60 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 61c 62 integer rowp(nshg,nnz), colm(nshg+1), 63 & iBC(nshg), 64 & ilwork(nlwork), 65 & iper(nshg), ifuncs(6) 66 67 real*8 vbc_prof(nshg,3) 68 69 integer stopjob 70 character*10 cname2 71 character*5 cname 72c 73c stuff for dynamic model s.w.avg and wall model 74c 75 dimension ifath(numnp), velbar(nfath,ndof), nsons(nfath) 76 77 dimension wallubar(2),walltot(2) 78c 79c.... For Farzin's Library 80c 81 integer eqnType, prjFlag, presPrjFlag, verbose 82c 83 real*8, allocatable, dimension(:,:) :: aperm, atemp, atempS 84 real*8, allocatable, dimension(:,:,:) :: apermS 85 86 real*8, allocatable, dimension(:,:) :: lhsP, lhsK, lhsS 87 real*8 almit, alfit, gamit 88c 89 character*1024 servername 90 character*20 fname1,fmt1 91 character*20 fname2,fmt2 92 character*60 fnamepold, fvarts 93 character*4 fname4c ! 4 characters 94 integer iarray(50) ! integers for headers 95 integer isgn(ndof), isgng(ndof) 96 97!MR CHANGE 98! real*8 rerr(nshg,10), ybar(nshg,13) ! including 7 sq. terms with 3 cross terms of uv, uw and vw 99! real*8 rerr(nshg,10), ybar(nshg,12) ! including 7 sq. terms with 3 cross terms of uv, uw and vw 100 real*8 rerr(nshg,10) 101 real*8, allocatable, dimension(:,:) :: ybar, strain, vorticity 102 real*8, allocatable, dimension(:,:) :: wallssVec, wallssVecbar 103!MR CHANGE 104 105 real*8 tcorecp(2), tcorecpscal(2) 106 107 integer, allocatable, dimension(:) :: ivarts 108 integer, allocatable, dimension(:) :: ivartsg 109 integer, allocatable, dimension(:) :: iv_rank 110 real*8, allocatable, dimension(:) :: vartssoln 111 real*8, allocatable, dimension(:) :: vartssolng 112 real*8, allocatable, dimension(:,:,:) :: yphbar 113 real*8 CFLworst(numel) 114 115!MR CHANGE 116 integer :: iv_rankpernode, iv_totnodes, iv_totcores 117 integer :: iv_node, iv_core, iv_thread 118!MR CHANGE 119 120 impistat = 0 121 impistat2 = 0 122 iISend = 0 123 iISendScal = 0 124 iIRecv = 0 125 iIRecvScal = 0 126 iWaitAll = 0 127 iWaitAllScal = 0 128 iAllR = 0 129 iAllRScal = 0 130 rISend = zero 131 rISendScal = zero 132 rIRecv = zero 133 rIRecvScal = zero 134 rWaitAll = zero 135 rWaitAllScal = zero 136 rAllR = zero 137 rAllRScal = zero 138 rCommu = zero 139 rCommuScal = zero 140 141 call initmemstat() 142 143c 144c hack SA variable 145c 146cHack BC(:,7)=BC(:,7)*0.001 147cHack if(lstep.eq.0) y(:,6)=y(:,6)*0.001 148 call SolverLicenseServer(servername) 149c 150c only master should be verbose 151c 152 153 if(numpe.gt.0 .and. myrank.ne.master)iverbose=0 154c 155 156 lskeep=lstep 157 158 inquire(file='xyzts.dat',exist=exts) 159 if(exts) then 160 161 open(unit=626,file='xyzts.dat',status='old') 162 read(626,*) ntspts, freq, tolpt, iterat, varcod 163 call sTD ! sets data structures 164 165 do jj=1,ntspts ! read coordinate data where solution desired 166 read(626,*) ptts(jj,1),ptts(jj,2),ptts(jj,3) 167 enddo 168 close(626) 169 170 statptts(:,:) = 0 171 parptts(:,:) = zero 172 varts(:,:) = zero 173 174 allocate (ivarts(ntspts*ndof)) 175 allocate (ivartsg(ntspts*ndof)) 176 allocate (iv_rank(ntspts)) 177 allocate (vartssoln(ntspts*ndof)) 178 allocate (vartssolng(ntspts*ndof)) 179 180 iv_rankpernode = iv_rankpercore*iv_corepernode 181 iv_totnodes = numpe/iv_rankpernode 182 iv_totcores = iv_corepernode*iv_totnodes 183 if (myrank .eq. 0) then 184 write(*,*) 'Info for probes:' 185 write(*,*) ' Ranks per core:',iv_rankpercore 186 write(*,*) ' Cores per node:',iv_corepernode 187 write(*,*) ' Ranks per node:',iv_rankpernode 188 write(*,*) ' Total number of nodes:',iv_totnodes 189 write(*,*) ' Total number of cores',iv_totcores 190 endif 191 192! if (myrank .eq. numpe-1) then 193 do jj=1,ntspts 194 195 ! Compute the adequate rank which will take care of probe jj 196 jjm1 = jj-1 197 iv_node = (iv_totnodes-1)-mod(jjm1,iv_totnodes) 198 iv_core = (iv_corepernode-1) - mod((jjm1 - 199 & mod(jjm1,iv_totnodes))/iv_totnodes,iv_corepernode) 200 iv_thread = (iv_rankpercore-1) - mod((jjm1- 201 & (mod(jjm1,iv_totcores)))/iv_totcores,iv_rankpercore) 202 iv_rank(jj) = iv_node*iv_rankpernode 203 & + iv_core*iv_rankpercore 204 & + iv_thread 205 206 if(myrank == 0) then 207 write(*,*) ' Probe', jj, 'handled by rank', 208 & iv_rank(jj), ' on node', iv_node 209 endif 210 211 ! Verification just in case 212 if(iv_rank(jj) .lt.0 .or. iv_rank(jj) .ge. numpe) then 213 write(*,*) 'WARNING: iv_rank(',jj,') is ', iv_rank(jj), 214 & ' and reset to numpe-1' 215 iv_rank(jj) = numpe-1 216 endif 217 218 ! Open the varts files 219 if(myrank == iv_rank(jj)) then 220 fvarts='varts/varts' 221 fvarts=trim(fvarts)//trim(cname2(jj)) 222 fvarts=trim(fvarts)//trim(cname2(lstep)) 223 fvarts=trim(fvarts)//'.dat' 224 fvarts=trim(fvarts) 225 open(unit=1000+jj, file=fvarts, status='unknown') 226 endif 227 enddo 228! endif 229 230 endif 231c 232c.... open history and aerodynamic forces files 233c 234 if (myrank .eq. master) then 235! open (unit=ihist, file=fhist, status='unknown') 236 open (unit=iforce, file=fforce, status='unknown') 237 open (unit=76, file="fort.76", status='unknown') 238 if(numImpSrfs.gt.0 .or. numRCRSrfs.gt.0) then 239 fnamepold = 'pold' 240 fnamepold = trim(fnamepold)//trim(cname2(lstep)) 241 fnamepold = trim(fnamepold)//'.dat' 242 fnamepold = trim(fnamepold) 243 open (unit=8177, file=fnamepold, status='unknown') 244 endif 245 endif 246c 247c.... initialize 248c 249 ifuncs(:) = 0 ! func. evaluation counter 250 istep = 0 251 yold = y 252 acold = ac 253 254 rerr = zero 255 256 if(ierrcalc.eq.1 .or. ioybar.eq.1) then ! we need ybar for error too 257 if (ivort == 1) then 258 allocate(ybar(nshg,17)) ! more space for vorticity if requested 259 else 260 allocate(ybar(nshg,13)) 261 endif 262 ybar = zero ! Initialize ybar to zero, which is essential 263 endif 264 265 if(ivort == 1) then 266 allocate(strain(nshg,6)) 267 allocate(vorticity(nshg,5)) 268 endif 269 270 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 271 allocate(wallssVec(nshg,3)) 272 if (ioybar .eq. 1) then 273 allocate(wallssVecbar(nshg,3)) 274 wallssVecbar = zero ! Initialization important if mean wss computed 275 endif 276 endif 277 278! both nstepsincycle and nphasesincycle needs to be set 279 if(nstepsincycle.eq.0) nphasesincycle = 0 280 if(nphasesincycle.ne.0) then 281! & allocate(yphbar(nshg,5,nphasesincycle)) 282 if (ivort == 1) then 283 allocate(yphbar(nshg,15,nphasesincycle)) ! more space for vorticity 284 else 285 allocate(yphbar(nshg,11,nphasesincycle)) 286 endif 287 yphbar = zero 288 endif 289 290!MR CHANGE END 291 292 vbc_prof(:,1:3) = BC(:,3:5) 293 if(iramp.eq.1) then 294 call BCprofileInit(vbc_prof,x) 295 endif 296 297c 298c.... ---------------> initialize Farzin's Library <--------------- 299c 300c.... assign parameter values 301c 302 do i = 1, 100 303 numeqns(i) = i 304 enddo 305 nKvecs = Kspace 306 prjFlag = iprjFlag 307 presPrjFlag = ipresPrjFlag 308 verbose = iverbose 309c 310c.... determine how many scalar equations we are going to need to solve 311c 312 nsolt=mod(impl(1),2) ! 1 if solving temperature 313 nsclrsol=nsolt+nsclr ! total number of scalars solved At 314 ! some point we probably want to create 315 ! a map, considering stepseq(), to find 316 ! what is actually solved and only 317 ! dimension lhs to the appropriate 318 ! size. (see 1.6.1 and earlier for a 319 ! "failed" attempt at this). 320 321 322 nsolflow=mod(impl(1),100)/10 ! 1 if solving flow 323 324c 325c.... Now, call Farzin's lesNew routine to initialize 326c memory space 327c 328 call genadj(colm, rowp, icnt ) ! preprocess the adjacency list 329 330 nnz_tot=icnt ! this is exactly the number of non-zero blocks on 331 ! this proc 332 333 if (nsolflow.eq.1) then 334 lesId = numeqns(1) 335 eqnType = 1 336 nDofs = 4 337 call myfLesNew( lesId, 41994, 338 & eqnType, 339 & nDofs, minIters, maxIters, 340 & nKvecs, prjFlag, nPrjs, 341 & presPrjFlag, nPresPrjs, epstol(1), 342 & prestol, verbose, statsflow, 343 & nPermDims, nTmpDims, servername ) 344 345 allocate (aperm(nshg,nPermDims)) 346 allocate (atemp(nshg,nTmpDims)) 347 allocate (lhsP(4,nnz_tot)) 348 allocate (lhsK(9,nnz_tot)) 349 350 call readLesRestart( lesId, aperm, nshg, myrank, lstep, 351 & nPermDims ) 352 353 else 354 nPermDims = 0 355 nTempDims = 0 356 endif 357 358 359 if(nsclrsol.gt.0) then 360 do isolsc=1,nsclrsol 361 lesId = numeqns(isolsc+1) 362 eqnType = 2 363 nDofs = 1 364 presPrjflag = 0 365 nPresPrjs = 0 366 prjFlag = 1 367 indx=isolsc+2-nsolt ! complicated to keep epstol(2) for 368 ! temperature followed by scalars 369 call myfLesNew( lesId, 41994, 370 & eqnType, 371 & nDofs, minIters, maxIters, 372 & nKvecs, prjFlag, nPrjs, 373 & presPrjFlag, nPresPrjs, epstol(indx), 374 & prestol, verbose, statssclr, 375 & nPermDimsS, nTmpDimsS, servername ) 376 enddo 377c 378c Assume all scalars have the same size needs 379c 380 allocate (apermS(nshg,nPermDimsS,nsclrsol)) 381 allocate (atempS(nshg,nTmpDimsS)) !they can all share this 382 allocate (lhsS(nnz_tot,nsclrsol)) 383c 384c actually they could even share with atemp but leave that for later 385c 386 else 387 nPermDimsS = 0 388 nTmpDimsS = 0 389 endif 390c 391c... prepare lumped mass if needed 392c 393 if((flmpr.ne.0).or.(flmpl.ne.0)) call genlmass(x, shp,shgl) 394c 395c.... -----------------> End of initialization <----------------- 396c 397c.....open the necessary files to gather time series 398c 399 lstep0 = lstep+1 400 nsteprcr = nstep(1)+lstep 401c 402c.... loop through the time sequences 403c 404 405 406 do 3000 itsq = 1, ntseq 407 itseq = itsq 408 409CAD tcorecp1 = second(0) 410CAD tcorewc1 = second(-1) 411c 412c.... set up the time integration parameters 413c 414 nstp = nstep(itseq) 415 nitr = niter(itseq) 416 LCtime = loctim(itseq) 417 dtol(:)= deltol(itseq,:) 418 419 call itrSetup ( y, acold ) 420c 421c...initialize the coefficients for the impedance convolution, 422c which are functions of alphaf so need to do it after itrSetup 423 if(numImpSrfs.gt.zero) then 424 call calcImpConvCoef (numImpSrfs, ntimeptpT) 425 endif 426c 427c...initialize the initial condition P(0)-RQ(0)-Pd(0) for RCR BC 428c need ndsurf so should be after initNABI 429 if(numRCRSrfs.gt.zero) then 430 call calcRCRic(y,nsrflistRCR,numRCRSrfs) 431 endif 432c 433c find the last solve of the flow in the step sequence so that we will 434c know when we are at/near end of step 435c 436c ilast=0 437 nitr=0 ! count number of flow solves in a step (# of iterations) 438 do i=1,seqsize 439 if(stepseq(i).eq.0) nitr=nitr+1 440 enddo 441 442 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 443 tcorecp(:) = zero ! used in solfar.f (solflow) 444 tcorecpscal(:) = zero ! used in solfar.f (solflow) 445 if(myrank.eq.0) then 446 tcorecp1 = TMRC() 447 endif 448 449c 450c.... loop through the time steps 451c 452 istop=0 453 rmub=datmat(1,2,1) 454 if(rmutarget.gt.0) then 455 rmue=rmutarget 456 else 457 rmue=datmat(1,2,1) ! keep constant 458 endif 459 460 if(iramp.eq.1) then 461 call BCprofileScale(vbc_prof,BC,yold) ! fix the yold values to the reset BC 462 isclr=1 ! fix scalar 463 do isclr=1,nsclr 464 call itrBCSclr(yold,ac,iBC,BC,iper,ilwork) 465 enddo 466 endif 467 468!MR CHANGE 469! do irank=0,numpe-1 470! if(irank==myrank) then 471! write(*,*) 'NUMTASK - rank,numtask: ',myrank,ilwork(1) 472! endif 473! call MPI_Barrier(MPI_COMM_WORLD,ierr) 474! enddo 475!MR CHANGE END 476 477 do 2000 istp = 1, nstp 478 if(iramp.eq.1) 479 & call BCprofileScale(vbc_prof,BC,yold) 480 481 call rerun_check(stopjob) 482 if(stopjob.ne.0) goto 2001 483 484 xi=istp*1.0/nstp 485 datmat(1,2,1)=rmub*(1.0-xi)+xi*rmue 486c write(*,*) "current mol. visc = ", datmat(1,2,1) 487c.... if we have time varying boundary conditions update the values of BC. 488c these will be for time step n+1 so use lstep+1 489c 490 if(itvn.gt.0) call BCint((lstep+1)*Delt(1), shp, shgl, 491 & shpb, shglb, x, BC, iBC) 492 493c 494c ... calculate the pressure contribution that depends on the history for the Imp. BC 495c 496 if(numImpSrfs.gt.0) then 497 call pHist(poldImp,QHistImp,ImpConvCoef, 498 & ntimeptpT,numImpSrfs) 499 if (myrank.eq.master) 500 & write(8177,*) (poldImp(n), n=1,numImpSrfs) 501 endif 502c 503c ... calc the pressure contribution that depends on the history for the RCR BC 504c 505 if(numRCRSrfs.gt.0) then 506 call CalcHopRCR (Delt(itseq), lstep, numRCRSrfs) 507 call CalcRCRConvCoef(lstep,numRCRSrfs) 508 call pHist(poldRCR,QHistRCR,RCRConvCoef,nsteprcr, 509 & numRCRSrfs) 510 if (myrank.eq.master) 511 & write(8177,*) (poldRCR(n), n=1,numRCRSrfs) 512 endif 513c 514c Decay of scalars 515c 516 if(nsclr.gt.0 .and. tdecay.ne.1) then 517 yold(:,6:ndof)=y(:,6:ndof)*tdecay 518 BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*tdecay 519 endif 520 521 if(nosource.eq.1) BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*0.8 522 523 524 if(iLES.gt.0) then !complicated stuff has moved to 525 !routine below 526 call lesmodels(yold, acold, shgl, shp, 527 & iper, ilwork, rowp, colm, 528 & nsons, ifath, x, 529 & iBC, BC) 530 531 532 endif 533 534c.... set traction BCs for modeled walls 535c 536 if (itwmod.ne.0) then 537 call asbwmod(yold, acold, x, BC, iBC, 538 & iper, ilwork, ifath, velbar) 539 endif 540 541!MR CHANGE 542c 543c.... Determine whether the vorticity field needs to be computed for this time step or not 544c 545 icomputevort = 0 546 if (ivort == 1) then ! Print vorticity = True in solver.inp 547 ! We then compute the vorticity only if we 548 ! 1) we write an intermediate checkpoint 549 ! 2) we reach the last time step and write the last checkpoint 550 ! 3) we accumulate statistics in ybar for every time step 551 ! BEWARE: we need here lstep+1 and istep+1 because the lstep and 552 ! istep gets incremened after the flowsolve, further below 553 if (((irs .ge. 1) .and. (mod(lstep+1, ntout) .eq. 0)) .or. 554 & istep+1.eq.nstep(itseq) .or. ioybar == 1) then 555 icomputevort = 1 556 endif 557 endif 558 559! write(*,*) 'icomputevort: ',icomputevort, ' - istep: ', 560! & istep,' - nstep(itseq):',nstep(itseq),'- lstep:', 561! & lstep, '- ntout:', ntout 562!MR CHANGE 563 564c 565c.... -----------------------> predictor phase <----------------------- 566c 567 call itrPredict(yold, y, acold, ac , uold, u, iBC) 568 call itrBC (y, ac, iBC, BC, iper,ilwork) 569 570 if(nsolt.eq.1) then 571 isclr=0 572 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 573 endif 574 do isclr=1,nsclr 575 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 576 enddo 577 iter=0 578 ilss=0 ! this is a switch thrown on first solve of LS redistance 579 do istepc=1,seqsize 580 icode=stepseq(istepc) 581 if(mod(icode,5).eq.0) then ! this is a solve 582 isolve=icode/10 583 if(icode.eq.0) then ! flow solve (encoded as 0) 584c 585 iter = iter+1 586 ifuncs(1) = ifuncs(1) + 1 587c 588 Force(1) = zero 589 Force(2) = zero 590 Force(3) = zero 591 HFlux = zero 592 lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 593 594 call SolFlow(y, ac, u, 595 & yold, acold, uold, 596 & x, iBC, 597 & BC, res, 598 & nPermDims, nTmpDims, aperm, 599 & atemp, iper, 600 & ilwork, shp, shgl, 601 & shpb, shglb, rowp, 602 & colm, lhsK, lhsP, 603 & solinc, rerr, tcorecp, 604 & GradV) 605 606 else ! scalar type solve 607 if (icode.eq.5) then ! Solve for Temperature 608 ! (encoded as (nsclr+1)*10) 609 isclr=0 610 ifuncs(2) = ifuncs(2) + 1 611 j=1 612 else ! solve a scalar (encoded at isclr*10) 613 isclr=isolve 614 ifuncs(isclr+2) = ifuncs(isclr+2) + 1 615 j=isclr+nsolt 616 if((iLSet.eq.2).and.(ilss.eq.0) 617 & .and.(isclr.eq.2)) then 618 ilss=1 ! throw switch (once per step) 619 y(:,7)=y(:,6) ! redistance field initialized 620 ac(:,7) = zero 621 call itrBCSclr ( y, ac, iBC, BC, iper, 622 & ilwork) 623c 624c....store the flow alpha, gamma parameter values and assigm them the 625c....Backward Euler parameters to solve the second levelset scalar 626c 627 alfit=alfi 628 gamit=gami 629 almit=almi 630 Deltt=Delt(1) 631 Dtglt=Dtgl 632 alfi = 1 633 gami = 1 634 almi = 1 635c Delt(1)= Deltt ! Give a pseudo time step 636 Dtgl = one / Delt(1) 637 endif ! level set eq. 2 638 endif ! deciding between temperature and scalar 639 640 lhs = 1 - min(1,mod(ifuncs(isclr+2)-1, 641 & LHSupd(isclr+2))) 642 643 call SolSclr(y, ac, u, 644 & yold, acold, uold, 645 & x, iBC, 646 & BC, nPermDimsS,nTmpDimsS, 647 & apermS(1,1,j), atempS, iper, 648 & ilwork, shp, shgl, 649 & shpb, shglb, rowp, 650 & colm, lhsS(1,j), 651 & solinc(1,isclr+5), tcorecpscal) 652 653 654 endif ! end of scalar type solve 655 656 else ! this is an update (mod did not equal zero) 657 iupdate=icode/10 ! what to update 658 if(icode.eq.1) then !update flow 659 call itrCorrect ( y, ac, u, solinc, iBC) 660 call itrBC (y, ac, iBC, BC, iper, ilwork) 661 else ! update scalar 662 isclr=iupdate !unless 663 if(icode.eq.6) isclr=0 664 if(iRANS.lt.-100)then ! RANS 665 call itrCorrectSclrPos(y,ac,solinc(1,isclr+5)) 666 else 667 call itrCorrectSclr (y, ac, solinc(1,isclr+5)) 668 endif 669 if (ilset.eq.2 .and. isclr.eq.2) then 670 if (ivconstraint .eq. 1) then 671 call itrBCSclr ( y, ac, iBC, BC, iper, 672 & ilwork) 673c 674c ... applying the volume constraint on second level set scalar 675c 676 call solvecon (y, x, iBC, BC, 677 & iper, ilwork, shp, shgl) 678c 679 endif ! end of volume constraint calculations 680 endif ! end of redistance calculations 681c 682 call itrBCSclr ( y, ac, iBC, BC, iper, 683 & ilwork) 684 endif ! end of flow or scalar update 685 endif ! end of switch between solve or update 686 enddo ! loop over sequence in step 687c 688c 689c.... obtain the time average statistics 690c 691 if (ioform .eq. 2) then 692 693 call stsGetStats( y, yold, ac, acold, 694 & u, uold, x, 695 & shp, shgl, shpb, shglb, 696 & iBC, BC, iper, ilwork, 697 & rowp, colm, lhsK, lhsP ) 698 endif 699 700c 701c Find the solution at the end of the timestep and move it to old 702c 703c 704c ...First to reassign the parameters for the original time integrator scheme 705c 706 if((iLSet.eq.2).and.(ilss.eq.1)) then 707 alfi =alfit 708 gami =gamit 709 almi =almit 710 Delt(1)=Deltt 711 Dtgl =Dtglt 712 endif 713 call itrUpdate( yold, acold, uold, y, ac, u) 714 call itrBC (yold, acold, iBC, BC, iper,ilwork) 715 716 istep = istep + 1 717 lstep = lstep + 1 718c 719c .. Print memory consumption on BGQ 720c 721 call printmeminfo("itrdrv"//char(0)) 722 723c 724c .. Compute vorticity 725c 726 if ( icomputevort == 1) then 727 728 ! vorticity components and magnitude 729 vorticity(:,1) = GradV(:,8)-GradV(:,6) !omega_x 730 vorticity(:,2) = GradV(:,3)-GradV(:,7) !omega_y 731 vorticity(:,3) = GradV(:,4)-GradV(:,2) !omega_z 732 vorticity(:,4) = sqrt( vorticity(:,1)*vorticity(:,1) 733 & + vorticity(:,2)*vorticity(:,2) 734 & + vorticity(:,3)*vorticity(:,3) ) 735 ! Q 736 strain(:,1) = GradV(:,1) !S11 737 strain(:,2) = 0.5*(GradV(:,2)+GradV(:,4)) !S12 738 strain(:,3) = 0.5*(GradV(:,3)+GradV(:,7)) !S13 739 strain(:,4) = GradV(:,5) !S22 740 strain(:,5) = 0.5*(GradV(:,6)+GradV(:,8)) !S23 741 strain(:,6) = GradV(:,9) !S33 742 743 vorticity(:,5) = 0.25*( vorticity(:,4)*vorticity(:,4) !Q 744 & - 2.0*( strain(:,1)*strain(:,1) 745 & + 2* strain(:,2)*strain(:,2) 746 & + 2* strain(:,3)*strain(:,3) 747 & + strain(:,4)*strain(:,4) 748 & + 2* strain(:,5)*strain(:,5) 749 & + strain(:,6)*strain(:,6))) 750 751 endif 752 753c 754c .. write out the solution 755c 756!MR CHANGE 757 758! if (((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) .or. 759! & lstep.eq.nstep(itseq)) then 760 if (((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) .or. 761 & istep.eq.nstep(itseq)) then 762!MR CHANGE END 763 764! Call to restar() will open restart file in write mode (and not append mode) 765! that is needed as other fields are written in append mode 766 call restar ('out ', yold ,ac) 767 if(ideformwall == 1) then 768 call write_displ(myrank, lstep, nshg, 3, uold ) 769 endif 770 771 if(ivort == 1) then 772 call write_field(myrank,'a','vorticity',9,vorticity, 773 & 'd',nshg,5,lstep) 774 endif 775!MR CHANGE 776! if(ioybar.eq.1) then 777! call write_field(myrank,'a','ybar',4, 778! & ybar,'d',nshg,13,lstep) 779! call write_field(myrank,'a','ybar',4, 780! & ybar,'d',nshg,12,lstep) 781! endif 782 783 call printmeminfo("itrdrv after checkpoint"//char(0)) 784!MR CHANGE END 785 786 endif 787 788c 789c ... update the flow history for the impedance convolution, filter it and write it out 790c 791 if(numImpSrfs.gt.zero) then 792 call UpdHistConv(y,nsrflistImp,numImpSrfs) !uses Delt(1) 793 endif 794 795c 796c ... update the flow history for the RCR convolution 797c 798 if(numRCRSrfs.gt.zero) then 799 call UpdHistConv(y,nsrflistRCR,numRCRSrfs) !uses lstep 800 endif 801 802c 803c.... compute the consistent boundary flux 804c 805 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 806 call Bflux ( yold, acold, uold, x, 807 & shp, shgl, shpb, 808 & shglb, ilwork, iBC, 809 & BC, iper, wallssVec) 810 endif 811 812c... dump TIME SERIES 813 814 if (exts) then 815 if (mod(lstep-1,freq).eq.0) then 816 817 if (numpe > 1) then 818 do jj = 1, ntspts 819 vartssoln((jj-1)*ndof+1:jj*ndof)=varts(jj,:) 820 ivarts=zero 821 enddo 822 do k=1,ndof*ntspts 823 if(vartssoln(k).ne.zero) ivarts(k)=1 824 enddo 825 826! call MPI_REDUCE(vartssoln, vartssolng, ndof*ntspts, 827! & MPI_DOUBLE_PRECISION, MPI_SUM, master, 828! & MPI_COMM_WORLD, ierr) 829 830 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 831 call MPI_ALLREDUCE(vartssoln, vartssolng, 832 & ndof*ntspts, 833 & MPI_DOUBLE_PRECISION, MPI_SUM, 834 & MPI_COMM_WORLD, ierr) 835 836! call MPI_REDUCE(ivarts, ivartsg, ndof*ntspts, 837! & MPI_INTEGER, MPI_SUM, master, 838! & MPI_COMM_WORLD, ierr) 839 840 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 841 call MPI_ALLREDUCE(ivarts, ivartsg, ndof*ntspts, 842 & MPI_INTEGER, MPI_SUM, 843 & MPI_COMM_WORLD, ierr) 844 845! if (myrank.eq.zero) then 846 do jj = 1, ntspts 847 848 if(myrank .eq. iv_rank(jj)) then 849 ! No need to update all varts components, only the one treated by the expected rank 850 ! Note: keep varts as a vector, as multiple probes could be treated by one rank 851 indxvarts = (jj-1)*ndof 852 do k=1,ndof 853 if(ivartsg(indxvarts+k).ne.0) then ! none of the vartssoln(parts) were non zero 854 varts(jj,k)=vartssolng(indxvarts+k)/ 855 & ivartsg(indxvarts+k) 856 endif 857 enddo 858 endif !only if myrank eq iv_rank(jj) 859 enddo 860! endif !only on master 861 endif !only if numpe > 1 862 863! if (myrank.eq.zero) then 864 do jj = 1, ntspts 865 if(myrank .eq. iv_rank(jj)) then 866 ifile = 1000+jj 867 write(ifile,555) lstep, (varts(jj,k),k=1,ndof) !Beware of format 555 - check ndof 868c call flush(ifile) 869 if (((irs .ge. 1) .and. 870 & (mod(lstep, ntout) .eq. 0))) then 871 close(ifile) 872 fvarts='varts/varts' 873 fvarts=trim(fvarts)//trim(cname2(jj)) 874 fvarts=trim(fvarts)//trim(cname2(lskeep)) 875 fvarts=trim(fvarts)//'.dat' 876 fvarts=trim(fvarts) 877 open(unit=ifile, file=fvarts, 878 & position='append') 879 endif !only when dumping restart 880 endif 881 enddo 882! endif !only on master 883 884 varts(:,:) = zero ! reset the array for next step 885 886 887!555 format(i6,5(2x,E12.5e2)) 888555 format(i6,6(2x,E20.12e2)) !assuming ndof = 6 here 889 890 endif 891 endif 892 893c 894c.... update and the aerodynamic forces 895c 896 call forces ( yold, ilwork ) 897 898 if((irscale.ge.0).or.(itwmod.gt.0)) 899 & call getvel (yold, ilwork, iBC, 900 & nsons, ifath, velbar) 901 902 if((irscale.ge.0).and.(myrank.eq.master)) then 903 call genscale(yold, x, iper, 904 & iBC, ifath, velbar, 905 & nsons) 906 endif 907c 908c.... print out results. 909c 910 ntoutv=max(ntout,100) ! velb is not needed so often 911 if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then 912 if( (mod(lstep, ntoutv) .eq. 0) .and. 913 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 914 & ((nsonmax.eq.1).and.(iLES.gt.0)))) 915 & call rwvelb ('out ', velbar ,ifail) 916 endif 917c 918c.... end of the NSTEP and NTSEQ loops 919c 920 921 922c 923c.... -------------------> error calculation <----------------- 924c 925 if(ierrcalc.eq.1 .or. ioybar.eq.1) then 926c$$$c 927c$$$c compute average 928c$$$c 929c$$$ tfact=one/istep 930c$$$ ybar =tfact*yold + (one-tfact)*ybar 931 932c compute average 933c ybar(:,1:3) are average velocity components 934c ybar(:,4) is average pressure 935c ybar(:,5) is average speed 936c ybar(:,6:8) is average of sq. of each vel. component 937c ybar(:,9) is average of sq. of pressure 938c ybar(:,10:12) is average of cross vel. components : uv, uw and vw 939c averaging procedure justified only for identical time step sizes 940c ybar(:,13) is average of eddy viscosity 941c ybar(:,14:16) is average vorticity components 942c ybar(:,17) is average vorticity magnitude 943c istep is number of time step 944c 945 icollectybar = 0 946 if(nphasesincycle.eq.0 .or. 947 & istep.gt.ncycles_startphaseavg*nstepsincycle) then 948 icollectybar = 1 949 if((istep-1).eq.ncycles_startphaseavg*nstepsincycle) 950 & istepsinybar = 0 ! init. to zero in first cycle in avg. 951 endif 952 953 if(icollectybar.eq.1) then 954 istepsinybar = istepsinybar+1 955 tfact=one/istepsinybar 956 957 if(myrank.eq.master .and. nphasesincycle.ne.0 .and. 958 & mod((istep-1),nstepsincycle).eq.0) 959 & write(*,*)'nsamples in phase average:',istepsinybar 960 961c ybar to contain the averaged ((u,v,w),p)-fields 962c and speed average, i.e., sqrt(u^2+v^2+w^2) 963c and avg. of sq. terms including 964c u^2, v^2, w^2, p^2 and cross terms of uv, uw and vw 965 966 ybar(:,1) = tfact*yold(:,1) + (one-tfact)*ybar(:,1) 967 ybar(:,2) = tfact*yold(:,2) + (one-tfact)*ybar(:,2) 968 ybar(:,3) = tfact*yold(:,3) + (one-tfact)*ybar(:,3) 969 ybar(:,4) = tfact*yold(:,4) + (one-tfact)*ybar(:,4) 970 ybar(:,5) = tfact*sqrt(yold(:,1)**2+yold(:,2)**2+ 971 & yold(:,3)**2) + (one-tfact)*ybar(:,5) 972 ybar(:,6) = tfact*yold(:,1)**2 + 973 & (one-tfact)*ybar(:,6) 974 ybar(:,7) = tfact*yold(:,2)**2 + 975 & (one-tfact)*ybar(:,7) 976 ybar(:,8) = tfact*yold(:,3)**2 + 977 & (one-tfact)*ybar(:,8) 978 ybar(:,9) = tfact*yold(:,4)**2 + 979 & (one-tfact)*ybar(:,9) 980 ybar(:,10) = tfact*yold(:,1)*yold(:,2) + !uv 981 & (one-tfact)*ybar(:,10) 982 ybar(:,11) = tfact*yold(:,1)*yold(:,3) + !uw 983 & (one-tfact)*ybar(:,11) 984 ybar(:,12) = tfact*yold(:,2)*yold(:,3) + !vw 985 & (one-tfact)*ybar(:,12) 986!MR CHANGE 987 if(nsclr.gt.0) !nut 988 & ybar(:,13) = tfact*yold(:,6) + (one-tfact)*ybar(:,13) 989 990 if(ivort == 1) then !vorticity 991 ybar(:,14) = tfact*vorticity(:,1) + 992 & (one-tfact)*ybar(:,14) 993 ybar(:,15) = tfact*vorticity(:,2) + 994 & (one-tfact)*ybar(:,15) 995 ybar(:,16) = tfact*vorticity(:,3) + 996 & (one-tfact)*ybar(:,16) 997 ybar(:,17) = tfact*vorticity(:,4) + 998 & (one-tfact)*ybar(:,17) 999 endif 1000 1001 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1002 wallssVecBar(:,1) = tfact*wallssVec(:,1) 1003 & +(one-tfact)*wallssVecBar(:,1) 1004 wallssVecBar(:,2) = tfact*wallssVec(:,2) 1005 & +(one-tfact)*wallssVecBar(:,2) 1006 wallssVecBar(:,3) = tfact*wallssVec(:,3) 1007 & +(one-tfact)*wallssVecBar(:,3) 1008 endif 1009!MR CHANGE END 1010 endif 1011c 1012c compute phase average 1013c 1014 if(nphasesincycle.ne.0 .and. 1015 & istep.gt.ncycles_startphaseavg*nstepsincycle) then 1016 1017c beginning of cycle is considered as ncycles_startphaseavg*nstepsincycle+1 1018 if((istep-1).eq.ncycles_startphaseavg*nstepsincycle) 1019 & icyclesinavg = 0 ! init. to zero in first cycle in avg. 1020 1021 ! find number of steps between phases 1022 nstepsbtwphase = nstepsincycle/nphasesincycle ! integer value 1023 if(mod(istep-1,nstepsincycle).eq.0) then 1024 iphase = 1 ! init. to one in beginning of every cycle 1025 icyclesinavg = icyclesinavg + 1 1026 endif 1027 1028 icollectphase = 0 1029 istepincycle = mod(istep,nstepsincycle) 1030 if(istepincycle.eq.0) istepincycle=nstepsincycle 1031 if(istepincycle.eq.iphase*nstepsbtwphase) then 1032 icollectphase = 1 1033 iphase = iphase+1 ! use 'iphase-1' below 1034 endif 1035 1036 if(icollectphase.eq.1) then 1037 tfactphase = one/icyclesinavg 1038 1039 if(myrank.eq.master) then 1040 write(*,*) 'nsamples in phase ',iphase-1,': ', 1041 & icyclesinavg 1042 endif 1043 1044 yphbar(:,1,iphase-1) = tfactphase*yold(:,1) + 1045 & (one-tfactphase)*yphbar(:,1,iphase-1) 1046 yphbar(:,2,iphase-1) = tfactphase*yold(:,2) + 1047 & (one-tfactphase)*yphbar(:,2,iphase-1) 1048 yphbar(:,3,iphase-1) = tfactphase*yold(:,3) + 1049 & (one-tfactphase)*yphbar(:,3,iphase-1) 1050 yphbar(:,4,iphase-1) = tfactphase*yold(:,4) + 1051 & (one-tfactphase)*yphbar(:,4,iphase-1) 1052 yphbar(:,5,iphase-1) = tfactphase*sqrt(yold(:,1)**2 1053 & +yold(:,2)**2+yold(:,3)**2) + 1054 & (one-tfactphase)*yphbar(:,5,iphase-1) 1055!MR CHANGE 1056 yphbar(:,6,iphase-1) = 1057 & tfactphase*yold(:,1)*yold(:,1) 1058 & +(one-tfactphase)*yphbar(:,6,iphase-1) 1059 1060 yphbar(:,7,iphase-1) = 1061 & tfactphase*yold(:,1)*yold(:,2) 1062 & +(one-tfactphase)*yphbar(:,7,iphase-1) 1063 1064 yphbar(:,8,iphase-1) = 1065 & tfactphase*yold(:,1)*yold(:,3) 1066 & +(one-tfactphase)*yphbar(:,8,iphase-1) 1067 1068 yphbar(:,9,iphase-1) = 1069 & tfactphase*yold(:,2)*yold(:,2) 1070 & +(one-tfactphase)*yphbar(:,9,iphase-1) 1071 1072 yphbar(:,10,iphase-1) = 1073 & tfactphase*yold(:,2)*yold(:,3) 1074 & +(one-tfactphase)*yphbar(:,10,iphase-1) 1075 1076 yphbar(:,11,iphase-1) = 1077 & tfactphase*yold(:,3)*yold(:,3) 1078 & +(one-tfactphase)*yphbar(:,11,iphase-1) 1079 1080 if(ivort == 1) then 1081 yphbar(:,12,iphase-1) = 1082 & tfactphase*vorticity(:,1) 1083 & +(one-tfactphase)*yphbar(:,12,iphase-1) 1084 yphbar(:,13,iphase-1) = 1085 & tfactphase*vorticity(:,2) 1086 & +(one-tfactphase)*yphbar(:,13,iphase-1) 1087 yphbar(:,14,iphase-1) = 1088 & tfactphase*vorticity(:,3) 1089 & +(one-tfactphase)*yphbar(:,14,iphase-1) 1090 yphbar(:,15,iphase-1) = 1091 & tfactphase*vorticity(:,4) 1092 & +(one-tfactphase)*yphbar(:,15,iphase-1) 1093 endif 1094!MR CHANGE END 1095 endif 1096 endif 1097c 1098c compute rms 1099c 1100 if(icollectybar.eq.1) then 1101 rerr(:, 7)=rerr(:, 7)+(yold(:,1)-ybar(:,1))**2 1102 rerr(:, 8)=rerr(:, 8)+(yold(:,2)-ybar(:,2))**2 1103 rerr(:, 9)=rerr(:, 9)+(yold(:,3)-ybar(:,3))**2 1104 rerr(:,10)=rerr(:,10)+(yold(:,4)-ybar(:,4))**2 1105 endif 1106 endif 1107 1108 if(istop.eq.1000) exit ! stop when delta small (see rstatic) 1109 2000 continue 1110 2001 continue 1111 1112 1113CAD tcorecp2 = second(0) 1114CAD tcorewc2 = second(-1) 1115 1116CAD write(6,*) 'T(core) cpu-wallclock = ',tcorecp2-tcorecp1, 1117CAD & tcorewc2-tcorewc1 1118 1119 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1120 if(myrank.eq.0) then 1121 tcorecp2 = TMRC() 1122 write(6,*) 'T(core) cpu = ',tcorecp2-tcorecp1 1123 write(6,*) '(Elm. form.',tcorecp(1), 1124 & ',Lin. alg. sol.',tcorecp(2),')' 1125 write(6,*) '(Elm. form. Scal.',tcorecpscal(1), 1126 & ',Lin. alg. sol. Scal.',tcorecpscal(2),')' 1127 write(6,*) '' 1128 1129 endif 1130 1131 call print_system_stats(tcorecp, tcorecpscal) 1132 call print_mesh_stats() 1133 call print_mpi_stats() 1134 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1135! return 1136c call MPI_Finalize() 1137c call MPI_ABORT(MPI_COMM_WORLD, ierr) 1138 1139 3000 continue 1140 1141c 1142c.... ----------------------> Post Processing <---------------------- 1143c 1144c.... print out the last step 1145c 1146 if ((irs .ge. 1) .and. ((mod(lstep, ntout) .ne. 0) .or. 1147 & (nstp .eq. 0))) then 1148 if( 1149 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 1150 & ((nsonmax.eq.1).and.iLES.gt.0))) 1151 & call rwvelb ('out ', velbar ,ifail) 1152 endif 1153 1154 1155 lesId = numeqns(1) 1156 1157!MR CHANGE 1158 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1159 if(myrank.eq.0) then 1160 tcormr1 = TMRC() 1161 endif 1162!MR CHANGE END 1163 1164 call saveLesRestart( lesId, aperm , nshg, myrank, lstep, 1165 & nPermDims ) 1166 1167!MR CHANGE 1168 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1169 if(myrank.eq.0) then 1170 tcormr2 = TMRC() 1171 write(6,*) 'Time to call saveLesRestart for projection and'// 1172 & 'pressure projection vectors', tcormr2-tcormr1 1173 endif 1174!MR CHANGE END 1175 1176 1177 if(ierrcalc.eq.1) then 1178c 1179c.....smooth the error indicators 1180c 1181 do i=1,ierrsmooth 1182 call errsmooth( rerr, x, iper, ilwork, shp, shgl, iBC ) 1183 end do 1184 1185!MR CHANGE 1186 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1187 if(myrank.eq.0) then 1188 tcormr1 = TMRC() 1189 endif 1190!MR CHANGE END 1191 1192 call write_error(myrank, lstep, nshg, 10, rerr ) 1193 1194!MR CHANGE 1195 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1196 if(myrank.eq.0) then 1197 call writeTimingMessage( 1198 & c_char_"error fields to " // c_null_char, 1199 & output_mode, 1200 & TMRC()-tcormr1) 1201 endif 1202!MR CHANGE END 1203 1204 1205 endif 1206 1207 if(ioybar.eq.1) then 1208 1209!MR CHANGE 1210! call write_field(myrank,'a','ybar',4, 1211! & ybar,'d',nshg,12,lstep) 1212 if(ivort == 1) then 1213 call write_field(myrank,'a','ybar',4, 1214 & ybar,'d',nshg,17,lstep) 1215 else 1216 call write_field(myrank,'a','ybar',4, 1217 & ybar,'d',nshg,13,lstep) 1218 endif 1219 deallocate(ybar) 1220 1221 1222 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1223 call write_field(myrank,'a','wssbar',6, 1224 & wallssVecBar,'d',nshg,3,lstep) 1225 deallocate(wallssVecbar) 1226 endif 1227 1228!MR CHANGE END 1229 if(nphasesincycle .gt. 0) then 1230 1231 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1232 if(myrank.eq.0) then 1233 tcormr1 = TMRC() 1234 endif 1235 1236 do iphase=1,nphasesincycle 1237 1238! call write_phavg(myrank,'w','phase average',13,iphase, 1239! & yphbar(:,:,iphase),'d',nshg,5,lstep) 1240! ! ybar field is repeated in files for phase average 1241! call write_phavg(myrank,'a','ybar',4,iphase, 1242! & ybar(:,1:5),'d',nshg,5,lstep) 1243 if(ivort == 1) then 1244 call write_phavg2(myrank,'a','phase_average',13,iphase, 1245 & nphasesincycle,yphbar(:,:,iphase),'d',nshg,15,lstep) 1246 else 1247 call write_phavg2(myrank,'a','phase_average',13,iphase, 1248 & nphasesincycle,yphbar(:,:,iphase),'d',nshg,11,lstep) 1249 endif 1250 1251 end do 1252 1253 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1254 if(myrank.eq.0) then 1255 tcormr2 = TMRC() 1256 write(6,*) 'Time to write all phase avg to the disks = ', 1257 & tcormr2-tcormr1 1258 endif 1259 deallocate(yphbar) 1260 endif 1261 1262 endif 1263 1264 if(ivort == 1) then 1265 deallocate(strain,vorticity) 1266 endif 1267 1268 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1269 deallocate(wallssVec) 1270 endif 1271 1272!MR CHANGE END 1273 1274 if ( ( ihessian .eq. 1 ) .and. ( numpe < 2 ) )then 1275 uhess = zero 1276 gradu = zero 1277 tf = zero 1278 1279 do ku=1,nshg 1280c tf(ku,1) = x(ku,1)**2+2*x(ku,1)*x(ku,2) 1281 tf(ku,1) = x(ku,1)**3 1282 end do 1283 1284 call hessian( yold, x, shp, shgl, iBC, 1285 & shpb, shglb, iper, ilwork, uhess, gradu ) 1286 1287 call write_hessian( uhess, gradu, nshg ) 1288 endif 1289 1290c if(iRANS.lt.0 .and. idistcalc.eq.1) then 1291 if(iRANS.lt.0) then 1292c call write_field(myrank,'a','DESd',4, 1293c & elDw,'d',numel,1,lstep) 1294 1295!MR CHANGE 1296 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1297 if(myrank.eq.0) then 1298 tcormr1 = TMRC() 1299 endif 1300!MR CHANGE END 1301 1302 call write_field(myrank,'a','dwal',4,d2wall,'d',nshg,1,lstep) 1303 deallocate(d2wall) 1304 1305!MR CHANGE 1306 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1307 if(myrank.eq.0) then 1308 tcormr2 = TMRC() 1309 if ( output_mode .eq. -1 ) then 1310 write(6,*) 'Time to write dwal to stream = ', tcormr2-tcormr1 1311 else 1312 write(6,*) 'Time to write dwal to disk = ', tcormr2-tcormr1 1313 endif 1314 endif 1315!MR CHANGE END 1316 1317 1318 endif 1319 1320c 1321c.... close history and aerodynamic forces files 1322c 1323 if (myrank .eq. master) then 1324! close (ihist) 1325 close (iforce) 1326 close(76) 1327 if(numImpSrfs.gt.0 .or. numRCRSrfs.gt.0) then 1328 close (8177) 1329 endif 1330 endif 1331c 1332c.... close varts file for probes 1333c 1334 if(exts) then 1335 do jj=1,ntspts 1336 if (myrank == iv_rank(jj)) then 1337 close(1000+jj) 1338 endif 1339 enddo 1340 deallocate (ivarts) 1341 deallocate (ivartsg) 1342 deallocate (iv_rank) 1343 deallocate (vartssoln) 1344 deallocate (vartssolng) 1345 endif 1346 1347!MR CHANGE 1348 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1349 if(myrank.eq.0) then 1350 write(*,*) 'itrdrv - done with aerodynamic forces' 1351 endif 1352!MR CHANGE 1353 1354 do isrf = 0,MAXSURF 1355! if ( nsrflist(isrf).ne.0 ) then 1356 if ( nsrflist(isrf).ne.0 .and. 1357 & myrank.eq.irankfilesforce(isrf)) then 1358 iunit=60+isrf 1359 close(iunit) 1360 endif 1361 enddo 1362 1363!MR CHANGE 1364 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1365 if(myrank.eq.0) then 1366 write(*,*) 'itrdrv - done with MAXSURF' 1367 endif 1368!MR CHANGE 1369 1370 1371 5 format(1X,F15.10,3X,F15.10,3X,F15.10,3X,F15.10) 1372 444 format(6(2x,e14.7)) 1373c 1374c.... end 1375c 1376 if(nsolflow.eq.1) then 1377 deallocate (lhsK) 1378 deallocate (lhsP) 1379 deallocate (aperm) 1380 deallocate (atemp) 1381 endif 1382 if(nsclrsol.gt.0) then 1383 deallocate (lhsS) 1384 deallocate (apermS) 1385 deallocate (atempS) 1386 endif 1387 1388 if(iabc==1) deallocate(acs) 1389 1390!MR CHANGE 1391 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1392 if(myrank.eq.0) then 1393 write(*,*) 'itrdrv - done - BACK TO process.f' 1394 endif 1395!MR CHANGE 1396 1397 1398 1399 return 1400 end 1401 1402 subroutine lesmodels(y, ac, shgl, shp, 1403 & iper, ilwork, rowp, colm, 1404 & nsons, ifath, x, 1405 & iBC, BC) 1406 1407 include "common.h" 1408 1409 real*8 y(nshg,ndof), ac(nshg,ndof), 1410 & x(numnp,nsd), 1411 & BC(nshg,ndofBC) 1412 real*8 shp(MAXTOP,maxsh,MAXQPT), 1413 & shgl(MAXTOP,nsd,maxsh,MAXQPT) 1414 1415c 1416 integer rowp(nshg,nnz), colm(nshg+1), 1417 & iBC(nshg), 1418 & ilwork(nlwork), 1419 & iper(nshg) 1420 dimension ifath(numnp), nsons(nfath) 1421 1422 real*8, allocatable, dimension(:) :: fwr2,fwr3,fwr4 1423 real*8, allocatable, dimension(:) :: stabdis,cdelsq1 1424 real*8, allocatable, dimension(:,:) :: xavegt, xavegt2,xavegt3 1425 1426 if( (iLES.gt.1) ) then ! Allocate Stuff for advanced LES models 1427 allocate (fwr2(nshg)) 1428 allocate (fwr3(nshg)) 1429 allocate (fwr4(nshg)) 1430 allocate (xavegt(nfath,12)) 1431 allocate (xavegt2(nfath,12)) 1432 allocate (xavegt3(nfath,12)) 1433 allocate (stabdis(nfath)) 1434 endif 1435 1436c.... get dynamic model coefficient 1437c 1438 ilesmod=iLES/10 1439c 1440c digit bit set filter rule, 10 bit set model 1441c 1442 if (ilesmod.eq.0) then ! 0 < iLES< 10 => dyn. model calculated 1443 ! at nodes based on discrete filtering 1444 1445 1446 if(isubmod.eq.2) then 1447 call SUPGdis(y, ac, shgl, shp, 1448 & iper, ilwork, 1449 & nsons, ifath, x, 1450 & iBC, BC, stabdis, xavegt3) 1451 endif 1452 1453 if( ((isubmod.eq.0).or.(isubmod.eq.2)))then ! If no 1454 ! sub-model 1455 ! or SUPG 1456 ! model wanted 1457 1458 if(i2filt.eq.0)then ! If simple filter 1459 1460 if(modlstats .eq. 0) then ! If no model stats wanted 1461 call getdmc (y, shgl, shp, 1462 & iper, ilwork, nsons, 1463 & ifath, x) 1464 else ! else get model stats 1465 call stdfdmc (y, shgl, shp, 1466 & iper, ilwork, nsons, 1467 & ifath, x) 1468 endif ! end of stats if statement 1469 1470 else ! else if twice filtering 1471 1472 call widefdmc(y, shgl, shp, 1473 & iper, ilwork, nsons, 1474 & ifath, x) 1475 1476 1477 endif ! end of simple filter if statement 1478 1479 endif ! end of SUPG or no sub-model if statement 1480 1481 1482 if( (isubmod.eq.1) ) then ! If DFWR sub-model wanted 1483 call cdelBHsq (y, shgl, shp, 1484 & iper, ilwork, nsons, 1485 & ifath, x, cdelsq1) 1486 call FiltRat (y, shgl, shp, 1487 & iper, ilwork, nsons, 1488 & ifath, x, cdelsq1, 1489 & fwr4, fwr3) 1490 1491 1492 if (i2filt.eq.0) then ! If simple filter wanted 1493 call DFWRsfdmc(y, shgl, shp, 1494 & iper, ilwork, nsons, 1495 & ifath, x, fwr2, fwr3) 1496 else ! else if twice filtering wanted 1497 call DFWRwfdmc(y, shgl, shp, 1498 & iper, ilwork, nsons, 1499 & ifath, x, fwr4, fwr4) 1500 endif ! end of simple filter if statement 1501 1502 endif ! end of DFWR sub-model if statement 1503 1504 if( (isubmod.eq.2) )then ! If SUPG sub-model wanted 1505 call dmcSUPG (y, ac, shgl, 1506 & shp, iper, ilwork, 1507 & nsons, ifath, x, 1508 & iBC, BC, rowp, colm, 1509 & xavegt2, stabdis) 1510 endif 1511 1512 if(idis.eq.1)then ! If SUPG/Model dissipation wanted 1513 call ediss (y, ac, shgl, 1514 & shp, iper, ilwork, 1515 & nsons, ifath, x, 1516 & iBC, BC, xavegt) 1517 endif 1518 1519 endif ! end of ilesmod 1520 1521 if (ilesmod .eq. 1) then ! 10 < iLES < 20 => dynamic-mixed 1522 ! at nodes based on discrete filtering 1523 call bardmc (y, shgl, shp, 1524 & iper, ilwork, 1525 & nsons, ifath, x) 1526 endif 1527 1528 if (ilesmod .eq. 2) then ! 20 < iLES < 30 => dynamic at quad 1529 ! pts based on lumped projection filt. 1530 1531 if(isubmod.eq.0)then 1532 call projdmc (y, shgl, shp, 1533 & iper, ilwork, x) 1534 else 1535 call cpjdmcnoi (y, shgl, shp, 1536 & iper, ilwork, x, 1537 & rowp, colm, 1538 & iBC, BC) 1539 endif 1540 1541 endif 1542 1543 if( (iLES.gt.1) ) then ! Deallocate Stuff for advanced LES models 1544 deallocate (fwr2) 1545 deallocate (fwr3) 1546 deallocate (fwr4) 1547 deallocate (xavegt) 1548 deallocate (xavegt2) 1549 deallocate (xavegt3) 1550 deallocate (stabdis) 1551 endif 1552 return 1553 end 1554 1555c 1556c...initialize the coefficients for the impedance convolution 1557c 1558 subroutine CalcImpConvCoef (numISrfs, numTpoints) 1559 1560 use convolImpFlow !uses flow history and impedance for convolution 1561 1562 include "common.h" !for alfi 1563 1564 integer numISrfs, numTpoints 1565 1566 allocate (ConvCoef(numTpoints+2,3)) !same time discret. for all imp. BC 1567 do j=1,numTpoints+2 1568 ConvCoef(j,:)=0.5/numTpoints !dt/2 divided by period T=N*dt 1569 ConvCoef(j,1)=ConvCoef(j,1)*(1.0-alfi)*(1.0-alfi) 1570 ConvCoef(j,2)=ConvCoef(j,2)*(1.0+2*alfi*(1.0-alfi)) 1571 ConvCoef(j,3)=ConvCoef(j,3)*alfi*alfi 1572 enddo 1573 ConvCoef(1,2)=zero 1574 ConvCoef(1,3)=zero 1575 ConvCoef(2,3)=zero 1576 ConvCoef(numTpoints+1,1)=zero 1577 ConvCoef(numTpoints+2,2)=zero 1578 ConvCoef(numTpoints+2,1)=zero 1579c 1580c...calculate the coefficients for the impedance convolution 1581c 1582 allocate (ImpConvCoef(numTpoints+2,numISrfs)) 1583 1584c..coefficients below assume Q linear in time step, Z constant 1585c do j=3,numTpoints 1586c ImpConvCoef(j,:) = ValueListImp(j-1,:)*ConvCoef(j,3) 1587c & + ValueListImp(j,:)*ConvCoef(j,2) 1588c & + ValueListImp(j+1,:)*ConvCoef(j,1) 1589c enddo 1590c ImpConvCoef(1,:) = ValueListImp(2,:)*ConvCoef(1,1) 1591c ImpConvCoef(2,:) = ValueListImp(2,:)*ConvCoef(2,2) 1592c & + ValueListImp(3,:)*ConvCoef(2,1) 1593c ImpConvCoef(numTpoints+1,:) = 1594c & ValueListImp(numTpoints,:)*ConvCoef(numTpoints+1,3) 1595c & + ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+1,2) 1596c ImpConvCoef(numTpoints+2,:) = 1597c & ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+2,3) 1598 1599c..try easiest convolution Q and Z constant per time step 1600 do j=3,numTpoints+1 1601 ImpConvCoef(j,:) = ValueListImp(j-1,:)/numTpoints 1602 enddo 1603 ImpConvCoef(1,:) =zero 1604 ImpConvCoef(2,:) =zero 1605 ImpConvCoef(numTpoints+2,:) = 1606 & ValueListImp(numTpoints+1,:)/numTpoints 1607c compensate for yalpha passed not y in Elmgmr() 1608 ImpConvCoef(numTpoints+1,:)= ImpConvCoef(numTpoints+1,:) 1609 & - ImpConvCoef(numTpoints+2,:)*(1.0-alfi)/alfi 1610 ImpConvCoef(numTpoints+2,:)= ImpConvCoef(numTpoints+2,:)/alfi 1611 return 1612 end 1613 1614c 1615c ... update the flow rate history for the impedance convolution, filter it and write it out 1616c 1617 subroutine UpdHistConv(y,nsrfIdList,numSrfs) 1618 1619 use convolImpFlow !brings ntimeptpT, QHistImp, QHistTry, QHistTryF, numImpSrfs 1620 use convolRCRFlow !brings QHistRCR, numRCRSrfs 1621 1622 include "common.h" 1623 1624 integer nsrfIdList(0:MAXSURF), numSrfs 1625 character*20 fname1 1626 character*10 cname2 1627 character*5 cname 1628 real*8 y(nshg,3) !velocity at time n+1 1629 real*8 NewQ(0:MAXSURF) !temporary unknown for the flow rate 1630 !that needs to be added to the flow history 1631 1632 call GetFlowQ(NewQ,y,nsrfIdList,numSrfs) !new flow at time n+1 1633c 1634c... for imp BC: shift QHist, add new constribution, filter and write out 1635c 1636 if(numImpSrfs.gt.zero) then 1637 do j=1, ntimeptpT 1638 QHistImp(j,1:numSrfs)=QHistImp(j+1,1:numSrfs) 1639 enddo 1640 QHistImp(ntimeptpT+1,1:numSrfs) = NewQ(1:numSrfs) 1641 1642c 1643c....filter the flow rate history 1644c 1645 cutfreq = 10 !hardcoded cutting frequency of the filter 1646 do j=1, ntimeptpT 1647 QHistTry(j,:)=QHistImp(j+1,:) 1648 enddo 1649 call Filter(QHistTryF,QHistTry,ntimeptpT,Delt(1),cutfreq) 1650c.... if no filter applied then uncomment next three lines 1651c do j=1, ntimeptpT 1652c QHistTryF(j,:)=QHistTry(j,:) 1653c enddo 1654 1655c QHistImp(1,:)=zero ! why do we do this? for beta(1,:) = zero it does not really matters 1656 do j=1, ntimeptpT 1657 QHistImp(j+1,:)=QHistTryF(j,:) 1658 enddo 1659c 1660c.... write out the new history of flow rates to Qhistor.dat 1661c 1662 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1663 & (istep .eq. nstep(1)))) .and. 1664 & (myrank .eq. master)) then 1665 open(unit=816, file='Qhistor.dat',status='replace') 1666 write(816,*) ntimeptpT 1667 do j=1,ntimeptpT+1 1668 write(816,*) (QHistImp(j,n),n=1, numSrfs) 1669 enddo 1670 close(816) 1671c... write out a copy with step number to be able to restart 1672 fname1 = 'Qhistor' 1673 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1674 open(unit=8166,file=trim(fname1),status='unknown') 1675 write(8166,*) ntimeptpT 1676 do j=1,ntimeptpT+1 1677 write(8166,*) (QHistImp(j,n),n=1, numSrfs) 1678 enddo 1679 close(8166) 1680 endif 1681 endif 1682 1683c 1684c... for RCR bc just add the new contribution 1685c 1686 if(numRCRSrfs.gt.zero) then 1687 QHistRCR(lstep+1,1:numSrfs) = NewQ(1:numSrfs) 1688c 1689c.... write out the new history of flow rates to Qhistor.dat 1690c 1691 if ((irs .ge. 1) .and. (myrank .eq. master)) then 1692 if(istep.eq.1) then 1693 open(unit=816,file='Qhistor.dat',status='unknown') 1694 else 1695 open(unit=816,file='Qhistor.dat',position='append') 1696 endif 1697 if(istep.eq.1) then 1698 do j=1,lstep 1699 write(816,*) j, (QHistRCR(j,n),n=1,numSrfs) ! read from file of previous run 1700 enddo 1701 endif 1702 write(816,*) lstep+1, (QHistRCR(lstep+1,n),n=1, numSrfs) 1703 close(816) 1704c... write out a copy with step number to be able to restart 1705 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1706 & (istep .eq. nstep(1)))) .and. 1707 & (myrank .eq. master)) then 1708 fname1 = 'Qhistor' 1709 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1710 open(unit=8166,file=trim(fname1),status='unknown') 1711 write(8166,*) lstep+1 1712 do j=1,lstep+1 1713 write(8166,*) (QHistRCR(j,n),n=1, numSrfs) 1714 enddo 1715 close(8166) 1716 endif 1717 endif 1718 endif 1719 1720 return 1721 end 1722 1723c 1724c...calculate the time varying coefficients for the RCR convolution 1725c 1726 subroutine CalcRCRConvCoef (stepn, numSrfs) 1727 1728 use convolRCRFlow !brings in ValueListRCR, dtRCR 1729 1730 include "common.h" !brings alfi 1731 1732 integer numSrfs, stepn 1733 1734 RCRConvCoef = zero 1735 if (stepn .eq. 0) then 1736 RCRConvCoef(1,:) = ValueListRCR(1,:)*(1.0-alfi) + 1737 & ValueListRCR(3,:)*(-alfi + 1.0 + 1/dtRCR(:) 1738 & - exp(-alfi*dtRCR(:))*(1 + 1/dtRCR(:))) 1739 RCRConvCoef(2,:) = ValueListRCR(1,:)*alfi 1740 & + ValueListRCR(3,:) 1741 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1742 endif 1743 if (stepn .ge. 1) then 1744 RCRConvCoef(1,:) =-ValueListRCR(3,:)*exp(-dtRCR(:)*(stepn+alfi)) 1745 & *(1 + (1 - exp(dtRCR(:)))/dtRCR(:)) 1746 RCRConvCoef(stepn+1,:) = ValueListRCR(1,:)*(1-alfi) 1747 & - ValueListRCR(3,:)*(alfi - 1 - 1/dtRCR(:) 1748 & + exp(-alfi*dtRCR(:))/dtRCR(:)*(2 - exp(-dtRCR(:)))) 1749 RCRConvCoef(stepn+2,:) = ValueListRCR(1,:)*alfi 1750 & + ValueListRCR(3,:) 1751 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1752 endif 1753 if (stepn .ge. 2) then 1754 do j=2,stepn 1755 RCRConvCoef(j,:) = ValueListRCR(3,:)/dtRCR(:)* 1756 & exp(-dtRCR(:)*(stepn + alfi + 2 - j))* 1757 & (1 - exp(dtRCR(:)))**2 1758 enddo 1759 endif 1760 1761c compensate for yalpha passed not y in Elmgmr() 1762 RCRConvCoef(stepn+1,:)= RCRConvCoef(stepn+1,:) 1763 & - RCRConvCoef(stepn+2,:)*(1.0-alfi)/alfi 1764 RCRConvCoef(stepn+2,:)= RCRConvCoef(stepn+2,:)/alfi 1765 1766 return 1767 end 1768 1769c 1770c...calculate the time dependent H operator for the RCR convolution 1771c 1772 subroutine CalcHopRCR (timestepRCR, stepn, numSrfs) 1773 1774 use convolRCRFlow !brings in HopRCR, dtRCR 1775 1776 include "common.h" 1777 1778 integer numSrfs, stepn 1779 real*8 PdistCur(0:MAXSURF), timestepRCR 1780 1781 HopRCR=zero 1782 call RCRint(timestepRCR*(stepn + alfi),PdistCur) 1783 HopRCR(1:numSrfs) = RCRic(1:numSrfs) 1784 & *exp(-dtRCR(1:numSrfs)*(stepn + alfi)) + PdistCur(1:numSrfs) 1785 return 1786 end 1787c 1788c ... initialize the influence of the initial conditions for the RCR BC 1789c 1790 subroutine calcRCRic(y,srfIdList,numSrfs) 1791 1792 use convolRCRFlow !brings RCRic, ValueListRCR, ValuePdist 1793 1794 include "common.h" 1795 1796 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled 1797 real*8 y(nshg,4) !need velocity and pressure 1798 real*8 Qini(0:MAXSURF) !initial flow rate 1799 real*8 PdistIni(0:MAXSURF) !initial distal pressure 1800 real*8 Pini(0:MAXSURF),CoupleArea(0:MAXSURF) ! initial pressure 1801 real*8 VelOnly(nshg,3), POnly(nshg) 1802 1803 allocate (RCRic(0:MAXSURF)) 1804 1805 if(lstep.eq.0) then 1806 VelOnly(:,1:3)=y(:,1:3) 1807 call GetFlowQ(Qini,VelOnly,srfIdList,numSrfs) !get initial flow 1808 QHistRCR(1,1:numSrfs)=Qini(1:numSrfs) !initialize QHistRCR 1809 POnly(:)=y(:,4) ! pressure 1810 call integrScalar(Pini,POnly,srfIdList,numSrfs) !get initial pressure integral 1811 POnly(:)=one ! one to get area 1812 call integrScalar(CoupleArea,POnly,srfIdList,numSrfs) !get surf area 1813 Pini(1:numSrfs) = Pini(1:numSrfs)/CoupleArea(1:numSrfs) 1814 else 1815 Qini(1:numSrfs)=QHistRCR(1,1:numSrfs) 1816 Pini(1:numSrfs)=zero ! hack 1817 endif 1818 call RCRint(istep,PdistIni) !get initial distal P (use istep) 1819 RCRic(1:numSrfs) = Pini(1:numSrfs) 1820 & - ValueListRCR(1,:)*Qini(1:numSrfs)-PdistIni(1:numSrfs) 1821 return 1822 end 1823 1824c.........function that integrates a scalar over a boundary 1825 subroutine integrScalar(scalInt,scal,srfIdList,numSrfs) 1826 1827 use pvsQbi !brings ndsurf, NASC 1828 1829 include "common.h" 1830 include "mpif.h" 1831 1832 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled, i, k 1833 real*8 scal(nshg), scalInt(0:MAXSURF), scalIntProc(0:MAXSURF) 1834 1835 scalIntProc = zero 1836 do i = 1,nshg 1837 if(numSrfs.gt.zero) then 1838 do k = 1,numSrfs 1839 irankCoupled = 0 1840 if (srfIdList(k).eq.ndsurf(i)) then 1841 irankCoupled=k 1842 scalIntProc(irankCoupled) = scalIntProc(irankCoupled) 1843 & + NASC(i)*scal(i) 1844 exit 1845 endif 1846 enddo 1847 endif 1848 enddo 1849c 1850c at this point, each scalint has its "nodes" contributions to the scalar 1851c accumulated into scalIntProc. Note, because NASC is on processor this 1852c will NOT be the scalar for the surface yet 1853c 1854c.... reduce integrated scalar for each surface, push on scalInt 1855c 1856 npars=MAXSURF+1 1857 call MPI_ALLREDUCE (scalIntProc, scalInt(:), npars, 1858 & MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr) 1859 1860 return 1861 end 1862 1863 subroutine writeTimingMessage(key,iomode,timing) 1864 use iso_c_binding 1865 use phstr 1866 implicit none 1867 1868 character(len=*) :: key 1869 integer :: iomode 1870 real*8 :: timing 1871 character(len=1024) :: 1872 & timing_msg = c_char_"Time to write "//c_null_char 1873 character(len=*), parameter :: 1874 & streamModeString = c_char_"stream"//c_null_char, 1875 & fileModeString = c_char_"disk"//c_null_char 1876 1877 call phstr_appendStr(timing_msg,key) 1878 if ( iomode .eq. -1 ) then 1879 call phstr_appendStr(timing_msg, streamModeString) 1880 else 1881 call phstr_appendStr(timing_msg, fileModeString) 1882 endif 1883 call phstr_appendStr(timing_msg, c_char_' = ' // c_null_char) 1884 call phstr_appendDbl(timing_msg, timing) 1885 write(6,*) timing_msg 1886 return 1887 end subroutine 1888 1889