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