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