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