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