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