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