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 use turbsa ! used to access d2wall 35 use wallData 36 use fncorpmod 37 use iso_c_binding 38 39c use readarrays !reads in uold and acold 40 41 include "common.h" 42 include "mpif.h" 43 include "auxmpi.h" 44#ifdef HAVE_SVLS 45 include "svLS.h" 46#endif 47#if !defined(HAVE_SVLS) && !defined(HAVE_LESLIB) 48#error "You must enable a linear solver during cmake setup" 49#endif 50 51c 52 53 54 real*8 y(nshg,ndof), ac(nshg,ndof), 55 & yold(nshg,ndof), acold(nshg,ndof), 56 & u(nshg,nsd), uold(nshg,nsd), 57 & x(numnp,nsd), solinc(nshg,ndof), 58 & BC(nshg,ndofBC), tf(nshg,ndof), 59 & GradV(nshg,nsdsq) 60 61c 62 real*8 res(nshg,ndof) 63c 64 real*8 shp(MAXTOP,maxsh,MAXQPT), 65 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 66 & shpb(MAXTOP,maxsh,MAXQPT), 67 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 68c 69 integer rowp(nshg,nnz), colm(nshg+1), 70 & iBC(nshg), 71 & ilwork(nlwork), 72 & iper(nshg), ifuncs(6) 73 74 real*8 vbc_prof(nshg,3) 75 76 integer stopjob 77 character*10 cname2 78 character*5 cname 79c 80c stuff for dynamic model s.w.avg and wall model 81c 82 dimension ifath(numnp), velbar(nfath,ndof), nsons(nfath) 83 84 dimension wallubar(2),walltot(2) 85c 86c.... For linear solver Library 87c 88 integer eqnType, prjFlag, presPrjFlag, verbose 89c 90 real*8, allocatable, dimension(:,:) :: aperm, atemp, atempS 91 real*8, allocatable, dimension(:,:,:) :: apermS 92 93 real*8, allocatable, dimension(:,:) :: lhsP, lhsK, lhsS 94 real*8 almit, alfit, gamit 95c 96 character*1024 servername 97 character*20 fname1,fmt1 98 character*20 fname2,fmt2 99 character*60 fnamepold, fvarts 100 character*4 fname4c ! 4 characters 101 integer iarray(50) ! integers for headers 102 integer isgn(ndof), isgng(ndof) 103 104!MR CHANGE 105! real*8 rerr(nshg,10), ybar(nshg,13) ! including 7 sq. terms with 3 cross terms of uv, uw and vw 106! real*8 rerr(nshg,10), ybar(nshg,12) ! including 7 sq. terms with 3 cross terms of uv, uw and vw 107 real*8 rerr(nshg,10) 108 real*8, allocatable, dimension(:,:) :: ybar, strain, vorticity 109 real*8, allocatable, dimension(:,:) :: wallssVec, wallssVecbar 110!MR CHANGE 111 112 real*8 tcorecp(2), tcorecpscal(2) 113 114 integer, allocatable, dimension(:) :: ivarts 115 integer, allocatable, dimension(:) :: ivartsg 116 integer, allocatable, dimension(:) :: iv_rank 117 real*8, allocatable, dimension(:) :: vartssoln 118 real*8, allocatable, dimension(:) :: vartssolng 119 real*8, allocatable, dimension(:,:,:) :: yphbar 120 real*8 CFLworst(numel) 121 122!MR CHANGE 123 integer :: iv_rankpernode, iv_totnodes, iv_totcores 124 integer :: iv_node, iv_core, iv_thread 125!MR CHANGE 126!-------------------------------------------------------------------- 127! Setting up svLS 128#ifdef HAVE_SVLS 129 INTEGER svLS_nFaces, gnNo, nNo, faIn, facenNo 130 INTEGER, ALLOCATABLE :: gNodes(:) 131 REAL*8, ALLOCATABLE :: sV(:,:) 132 133 TYPE(svLS_commuType) communicator 134 TYPE(svLS_lhsType) svLS_lhs 135 TYPE(svLS_lsType) svLS_ls 136! repeat for scalar solves (up to 4 at this time which is consistent with rest of PHASTA) 137 TYPE(svLS_commuType) communicator_S(4) 138 TYPE(svLS_lhsType) svLS_lhs_S(4) 139 TYPE(svLS_lsType) svLS_ls_S(4) 140#endif 141 142 impistat = 0 143 impistat2 = 0 144 iISend = 0 145 iISendScal = 0 146 iIRecv = 0 147 iIRecvScal = 0 148 iWaitAll = 0 149 iWaitAllScal = 0 150 iAllR = 0 151 iAllRScal = 0 152 rISend = zero 153 rISendScal = zero 154 rIRecv = zero 155 rIRecvScal = zero 156 rWaitAll = zero 157 rWaitAllScal = zero 158 rAllR = zero 159 rAllRScal = zero 160 rCommu = zero 161 rCommuScal = zero 162 163 call initmemstat() 164 165c 166c hack SA variable 167c 168cHack BC(:,7)=BC(:,7)*0.001 169cHack if(lstep.eq.0) y(:,6)=y(:,6)*0.001 170!-------------------------------------------------------------------- 171! Setting up svLS Moved down for better org 172 173#ifdef HAVE_LESLIB 174! IF (svLSFlag .EQ. 0) THEN !When we get a PETSc option it also could block this or a positive leslib 175 call SolverLicenseServer(servername) 176! ENDIF 177#endif 178c 179c only master should be verbose 180c 181 182 if(numpe.gt.0 .and. myrank.ne.master)iverbose=0 183c 184 185 lskeep=lstep 186 187 inquire(file='xyzts.dat',exist=exts) 188 if(exts) then 189 190 open(unit=626,file='xyzts.dat',status='old') 191 read(626,*) ntspts, freq, tolpt, iterat, varcod 192 call sTD ! sets data structures 193 194 do jj=1,ntspts ! read coordinate data where solution desired 195 read(626,*) ptts(jj,1),ptts(jj,2),ptts(jj,3) 196 enddo 197 close(626) 198 199 statptts(:,:) = 0 200 parptts(:,:) = zero 201 varts(:,:) = zero 202 203 allocate (ivarts(ntspts*ndof)) 204 allocate (ivartsg(ntspts*ndof)) 205 allocate (iv_rank(ntspts)) 206 allocate (vartssoln(ntspts*ndof)) 207 allocate (vartssolng(ntspts*ndof)) 208 209 iv_rankpernode = iv_rankpercore*iv_corepernode 210 iv_totnodes = numpe/iv_rankpernode 211 iv_totcores = iv_corepernode*iv_totnodes 212 if (myrank .eq. 0) then 213 write(*,*) 'Info for probes:' 214 write(*,*) ' Ranks per core:',iv_rankpercore 215 write(*,*) ' Cores per node:',iv_corepernode 216 write(*,*) ' Ranks per node:',iv_rankpernode 217 write(*,*) ' Total number of nodes:',iv_totnodes 218 write(*,*) ' Total number of cores',iv_totcores 219 endif 220 221! if (myrank .eq. numpe-1) then 222 do jj=1,ntspts 223 224 ! Compute the adequate rank which will take care of probe jj 225 jjm1 = jj-1 226 iv_node = (iv_totnodes-1)-mod(jjm1,iv_totnodes) 227 iv_core = (iv_corepernode-1) - mod((jjm1 - 228 & mod(jjm1,iv_totnodes))/iv_totnodes,iv_corepernode) 229 iv_thread = (iv_rankpercore-1) - mod((jjm1- 230 & (mod(jjm1,iv_totcores)))/iv_totcores,iv_rankpercore) 231 iv_rank(jj) = iv_node*iv_rankpernode 232 & + iv_core*iv_rankpercore 233 & + iv_thread 234 235 if(myrank == 0) then 236 write(*,*) ' Probe', jj, 'handled by rank', 237 & iv_rank(jj), ' on node', iv_node 238 endif 239 240 ! Verification just in case 241 if(iv_rank(jj) .lt.0 .or. iv_rank(jj) .ge. numpe) then 242 write(*,*) 'WARNING: iv_rank(',jj,') is ', iv_rank(jj), 243 & ' and reset to numpe-1' 244 iv_rank(jj) = numpe-1 245 endif 246 247 ! Open the varts files 248 if(myrank == iv_rank(jj)) then 249 fvarts='varts/varts' 250 fvarts=trim(fvarts)//trim(cname2(jj)) 251 fvarts=trim(fvarts)//trim(cname2(lstep)) 252 fvarts=trim(fvarts)//'.dat' 253 fvarts=trim(fvarts) 254 open(unit=1000+jj, file=fvarts, status='unknown') 255 endif 256 enddo 257! endif 258 259 endif 260c 261c.... open history and aerodynamic forces files 262c 263 if (myrank .eq. master) then 264 open (unit=ihist, file=fhist, status='unknown') 265 open (unit=iforce, file=fforce, status='unknown') 266 open (unit=76, file="fort.76", status='unknown') 267 if(numImpSrfs.gt.0 .or. numRCRSrfs.gt.0) then 268 fnamepold = 'pold' 269 fnamepold = trim(fnamepold)//trim(cname2(lstep)) 270 fnamepold = trim(fnamepold)//'.dat' 271 fnamepold = trim(fnamepold) 272 open (unit=8177, file=fnamepold, status='unknown') 273 endif 274 endif 275c 276c.... initialize 277c 278 ifuncs(:) = 0 ! func. evaluation counter 279 istep = 0 280 yold = y 281 acold = ac 282 283 rerr = zero 284 285 if(ierrcalc.eq.1 .or. ioybar.eq.1) then ! we need ybar for error too 286 if (ivort == 1) then 287 allocate(ybar(nshg,17)) ! more space for vorticity if requested 288 else 289 allocate(ybar(nshg,13)) 290 endif 291 ybar = zero ! Initialize ybar to zero, which is essential 292 endif 293 294 if(ivort == 1) then 295 allocate(strain(nshg,6)) 296 allocate(vorticity(nshg,5)) 297 endif 298 299 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 300 allocate(wallssVec(nshg,3)) 301 if (ioybar .eq. 1) then 302 allocate(wallssVecbar(nshg,3)) 303 wallssVecbar = zero ! Initialization important if mean wss computed 304 endif 305 endif 306 307! both nstepsincycle and nphasesincycle needs to be set 308 if(nstepsincycle.eq.0) nphasesincycle = 0 309 if(nphasesincycle.ne.0) then 310! & allocate(yphbar(nshg,5,nphasesincycle)) 311 if (ivort == 1) then 312 allocate(yphbar(nshg,15,nphasesincycle)) ! more space for vorticity 313 else 314 allocate(yphbar(nshg,11,nphasesincycle)) 315 endif 316 yphbar = zero 317 endif 318 319!MR CHANGE END 320 321 vbc_prof(:,1:3) = BC(:,3:5) 322 if(iramp.eq.1) then 323 call BCprofileInit(vbc_prof,x) 324 endif 325 326c 327c.... ---------------> initialize LesLib Library <--------------- 328c 329c.... assign parameter values 330c 331 do i = 1, 100 332 numeqns(i) = i 333 enddo 334 nKvecs = Kspace 335 prjFlag = iprjFlag 336 presPrjFlag = ipresPrjFlag 337 verbose = iverbose 338c 339c.... determine how many scalar equations we are going to need to solve 340c 341 nsolt=mod(impl(1),2) ! 1 if solving temperature 342 nsclrsol=nsolt+nsclr ! total number of scalars solved At 343 ! some point we probably want to create 344 ! a map, considering stepseq(), to find 345 ! what is actually solved and only 346 ! dimension lhs to the appropriate 347 ! size. (see 1.6.1 and earlier for a 348 ! "failed" attempt at this). 349 350 351 nsolflow=mod(impl(1),100)/10 ! 1 if solving flow 352 353c 354c.... Now, call lesNew routine to initialize 355c memory space 356c 357 call genadj(colm, rowp, icnt ) ! preprocess the adjacency list 358 359 nnz_tot=icnt ! this is exactly the number of non-zero blocks on 360 ! this proc 361 362 if (nsolflow.eq.1) then ! start of setup for the flow 363 lesId = numeqns(1) 364 eqnType = 1 365 nDofs = 4 366 367!-------------------------------------------------------------------- 368! Rest of configuration of svLS is added here, where we have LHS 369! pointers 370 371 nPermDims = 1 372 nTmpDims = 1 373 374 allocate (lhsP(4,nnz_tot)) 375 allocate (lhsK(9,nnz_tot)) 376 377! Setting up svLS or leslib for flow 378 379 IF (svLSFlag .EQ. 1) THEN 380#ifdef HAVE_SVLS 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#else 447 if(myrank.eq.0) write(*,*) 'your input requests svLS but your cmake did not build for it' 448 call error('itrdrv ','nosVLS',svLSFlag) 449#endif 450 ENDIF 451 452 if(leslib.eq.1) then 453#ifdef HAVE_LESLIB 454!-------------------------------------------------------------------- 455 call myfLesNew( lesId, 41994, 456 & eqnType, 457 & nDofs, minIters, maxIters, 458 & nKvecs, prjFlag, nPrjs, 459 & presPrjFlag, nPresPrjs, epstol(1), 460 & prestol, verbose, statsflow, 461 & nPermDims, nTmpDims, servername ) 462 463 allocate (aperm(nshg,nPermDims)) 464 allocate (atemp(nshg,nTmpDims)) 465 call readLesRestart( lesId, aperm, nshg, myrank, lstep, 466 & nPermDims ) 467#else 468 if(myrank.eq.0) write(*,*) 'your input requests leslib but your cmake did not build for it' 469 call error('itrdrv ','nolslb',leslib) 470#endif 471 endif !leslib=1 472 473 else ! not solving flow just scalar 474 nPermDims = 0 475 nTempDims = 0 476 endif 477 478 479 if(nsclrsol.gt.0) then 480 do isolsc=1,nsclrsol 481 lesId = numeqns(isolsc+1) 482 eqnType = 2 483 nDofs = 1 484 presPrjflag = 0 485 nPresPrjs = 0 486 prjFlag = 1 487 indx=isolsc+2-nsolt ! complicated to keep epstol(2) for 488 ! temperature followed by scalars 489! Setting up svLS or leslib for scalar 490#ifdef HAVE_SVLS 491 IF (svLSFlag .EQ. 1) THEN 492 svLSType=2 !only option for scalars 493! reltol for the GMRES is the stop criterion 494! also using Kspaceand maxIters from setup for ACUSIM 495! 496 CALL svLS_LS_CREATE(svLS_ls_S(isolsc), svLSType, dimKry=Kspace, 497 2 relTol=epstol(indx), 498 3 maxItr=maxIters 499 4 ) 500 501 CALL svLS_COMMU_CREATE(communicator_S(isolsc), MPI_COMM_WORLD) 502 503 svLS_nFaces = 1 !not sure about this...should try it with zero 504 505 CALL svLS_LHS_CREATE(svLS_lhs_S(isolsc), communicator_S(isolsc), gnNo, nNo, 506 2 nnz_tot, ltg, colm, rowp, svLS_nFaces) 507 508 faIn = 1 509 facenNo = 0 510 ib=5+isolsc 511 DO i=1, nshg 512 IF (btest(iBC(i),ib)) facenNo = facenNo + 1 513 END DO 514 ALLOCATE(gNodes(facenNo), sV(1,facenNo)) 515 sV = 0D0 516 j = 0 517 DO i=1, nshg 518 IF (btest(iBC(i),ib)) THEN 519 j = j + 1 520 gNodes(j) = i 521 END IF 522 END DO 523 524 CALL svLS_BC_CREATE(svLS_lhs_S(isolsc), faIn, facenNo, 525 2 1, BC_TYPE_Dir, gNodes, sV(1,:)) 526 DEALLOCATE(gNodes) 527 DEALLOCATE(sV) 528 529 if( isolsc.eq.1) then ! if multiple scalars make sure done once 530 allocate (apermS(1,1,1)) 531 allocate (atempS(1,1)) !they can all share this 532 endif 533 ENDIF !svLS handing scalar solve 534#endif 535 536 537#ifdef HAVE_LESLIB 538 if (leslib.eq.1) then 539 call myfLesNew( lesId, 41994, 540 & eqnType, 541 & nDofs, minIters, maxIters, 542 & nKvecs, prjFlag, nPrjs, 543 & presPrjFlag, nPresPrjs, epstol(indx), 544 & prestol, verbose, statssclr, 545 & nPermDimsS, nTmpDimsS, servername ) 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_LESLIB 551 if(leslib.eq.1) 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.... update and the aerodynamic forces 945c 946 call forces ( yold, ilwork ) 947 948c 949c .. write out the instantaneous solution 950c 9512001 continue ! we could get here by 2001 label if user requested stop 952 if (((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) .or. 953 & istep.eq.nstep(itseq)) then 954 955!so that we can see progress in force file close it so that it flushes 956!and then reopen in append mode 957 958 close(iforce) 959 open (unit=iforce, file=fforce, position='append') 960 961! Call to restar() will open restart file in write mode (and not append mode) 962! that is needed as other fields are written in append mode 963 964 call restar ('out ', yold ,ac) 965 if(ideformwall == 1) then 966! call write_displ(myrank, lstep, nshg, 3, uold ) 967 if(myrank.eq.master) then 968 write(*,*) 'ideformwall is 1 and is a dead code path... exiting' 969 endif 970 stop 971 endif 972 973 if(ivort == 1) then 974 call write_field(myrank,'a','vorticity',9,vorticity, 975 & 'd',nshg,5,lstep) 976 endif 977 978 call printmeminfo("itrdrv after checkpoint"//char(0)) 979 else if(stopjob.eq.-2) then 980 if(myrank.eq.master) then 981 write(*,*) 'line 755 says no write before stopping' 982 write(*,*) 'istep,nstep,irs',istep,nstep(itseq),irs 983 endif 984 endif !just the instantaneous stuff for videos 985c 986c.... compute the consistent boundary flux 987c 988 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 989 call Bflux ( yold, acold, uold, x, 990 & shp, shgl, shpb, 991 & shglb, ilwork, iBC, 992 & BC, iper, wallssVec) 993 endif 994 995 if(stopjob.eq.-2) goto 2003 996 997 998c 999c ... update the flow history for the impedance convolution, filter it and write it out 1000c 1001 if(numImpSrfs.gt.zero) then 1002 call UpdHistConv(y,nsrflistImp,numImpSrfs) !uses Delt(1) 1003 endif 1004 1005c 1006c ... update the flow history for the RCR convolution 1007c 1008 if(numRCRSrfs.gt.zero) then 1009 call UpdHistConv(y,nsrflistRCR,numRCRSrfs) !uses lstep 1010 endif 1011 1012 1013c... dump TIME SERIES 1014 1015 if (exts) then 1016 if (mod(lstep-1,freq).eq.0) then 1017 1018 if (numpe > 1) then 1019 do jj = 1, ntspts 1020 vartssoln((jj-1)*ndof+1:jj*ndof)=varts(jj,:) 1021 ivarts=zero 1022 enddo 1023 do k=1,ndof*ntspts 1024 if(vartssoln(k).ne.zero) ivarts(k)=1 1025 enddo 1026 1027! call MPI_REDUCE(vartssoln, vartssolng, ndof*ntspts, 1028! & MPI_DOUBLE_PRECISION, MPI_SUM, master, 1029! & MPI_COMM_WORLD, ierr) 1030 1031 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1032 call MPI_ALLREDUCE(vartssoln, vartssolng, 1033 & ndof*ntspts, 1034 & MPI_DOUBLE_PRECISION, MPI_SUM, 1035 & MPI_COMM_WORLD, ierr) 1036 1037! call MPI_REDUCE(ivarts, ivartsg, ndof*ntspts, 1038! & MPI_INTEGER, MPI_SUM, master, 1039! & MPI_COMM_WORLD, ierr) 1040 1041 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1042 call MPI_ALLREDUCE(ivarts, ivartsg, ndof*ntspts, 1043 & MPI_INTEGER, MPI_SUM, 1044 & MPI_COMM_WORLD, ierr) 1045 1046! if (myrank.eq.zero) then 1047 do jj = 1, ntspts 1048 1049 if(myrank .eq. iv_rank(jj)) then 1050 ! No need to update all varts components, only the one treated by the expected rank 1051 ! Note: keep varts as a vector, as multiple probes could be treated by one rank 1052 indxvarts = (jj-1)*ndof 1053 do k=1,ndof 1054 if(ivartsg(indxvarts+k).ne.0) then ! none of the vartssoln(parts) were non zero 1055 varts(jj,k)=vartssolng(indxvarts+k)/ 1056 & ivartsg(indxvarts+k) 1057 endif 1058 enddo 1059 endif !only if myrank eq iv_rank(jj) 1060 enddo 1061! endif !only on master 1062 endif !only if numpe > 1 1063 1064! if (myrank.eq.zero) then 1065 do jj = 1, ntspts 1066 if(myrank .eq. iv_rank(jj)) then 1067 ifile = 1000+jj 1068 write(ifile,555) lstep, (varts(jj,k),k=1,ndof) !Beware of format 555 - check ndof 1069c call flush(ifile) 1070 if (((irs .ge. 1) .and. 1071 & (mod(lstep, ntout) .eq. 0))) then 1072 close(ifile) 1073 fvarts='varts/varts' 1074 fvarts=trim(fvarts)//trim(cname2(jj)) 1075 fvarts=trim(fvarts)//trim(cname2(lskeep)) 1076 fvarts=trim(fvarts)//'.dat' 1077 fvarts=trim(fvarts) 1078 open(unit=ifile, file=fvarts, 1079 & position='append') 1080 endif !only when dumping restart 1081 endif 1082 enddo 1083! endif !only on master 1084 1085 varts(:,:) = zero ! reset the array for next step 1086 1087 1088!555 format(i6,5(2x,E12.5e2)) 1089555 format(i6,6(2x,E20.12e2)) !assuming ndof = 6 here 1090 1091 endif 1092 endif 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 call destroyWallData 1484 call destroyfncorp 1485 1486 3000 continue 1487 1488 1489c 1490c.... close history and aerodynamic forces files 1491c 1492 if (myrank .eq. master) then 1493! close (ihist) 1494 close (iforce) 1495 close(76) 1496 if(numImpSrfs.gt.0 .or. numRCRSrfs.gt.0) then 1497 close (8177) 1498 endif 1499 endif 1500c 1501c.... close varts file for probes 1502c 1503 if(exts) then 1504 do jj=1,ntspts 1505 if (myrank == iv_rank(jj)) then 1506 close(1000+jj) 1507 endif 1508 enddo 1509 deallocate (ivarts) 1510 deallocate (ivartsg) 1511 deallocate (iv_rank) 1512 deallocate (vartssoln) 1513 deallocate (vartssolng) 1514 endif 1515 1516!MR CHANGE 1517 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1518 if(myrank.eq.0) then 1519 write(*,*) 'itrdrv - done with aerodynamic forces' 1520 endif 1521!MR CHANGE 1522 1523 do isrf = 0,MAXSURF 1524! if ( nsrflist(isrf).ne.0 ) then 1525 if ( nsrflist(isrf).ne.0 .and. 1526 & myrank.eq.irankfilesforce(isrf)) then 1527 iunit=60+isrf 1528 close(iunit) 1529 endif 1530 enddo 1531 1532!MR CHANGE 1533 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1534 if(myrank.eq.0) then 1535 write(*,*) 'itrdrv - done with MAXSURF' 1536 endif 1537!MR CHANGE 1538 1539 1540 5 format(1X,F15.10,3X,F15.10,3X,F15.10,3X,F15.10) 1541 444 format(6(2x,e14.7)) 1542c 1543c.... end 1544c 1545 if(nsolflow.eq.1) then 1546 deallocate (lhsK) 1547 deallocate (lhsP) 1548 IF (svLSFlag .NE. 1) THEN 1549 deallocate (aperm) 1550 deallocate (atemp) 1551 ENDIF 1552 endif 1553 if(nsclrsol.gt.0) then 1554 deallocate (lhsS) 1555 deallocate (apermS) 1556 deallocate (atempS) 1557 endif 1558 1559 if(iabc==1) deallocate(acs) 1560 1561!MR CHANGE 1562 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1563 if(myrank.eq.0) then 1564 write(*,*) 'itrdrv - done - BACK TO process.f' 1565 endif 1566!MR CHANGE 1567 1568 1569 1570 return 1571 end 1572 1573 subroutine lesmodels(y, ac, shgl, shp, 1574 & iper, ilwork, rowp, colm, 1575 & nsons, ifath, x, 1576 & iBC, BC) 1577 1578 include "common.h" 1579 1580 real*8 y(nshg,ndof), ac(nshg,ndof), 1581 & x(numnp,nsd), 1582 & BC(nshg,ndofBC) 1583 real*8 shp(MAXTOP,maxsh,MAXQPT), 1584 & shgl(MAXTOP,nsd,maxsh,MAXQPT) 1585 1586c 1587 integer rowp(nshg,nnz), colm(nshg+1), 1588 & iBC(nshg), 1589 & ilwork(nlwork), 1590 & iper(nshg) 1591 dimension ifath(numnp), nsons(nfath) 1592 1593 real*8, allocatable, dimension(:) :: fwr2,fwr3,fwr4 1594 real*8, allocatable, dimension(:) :: stabdis,cdelsq1 1595 real*8, allocatable, dimension(:,:) :: xavegt, xavegt2,xavegt3 1596 1597 if( (iLES.gt.1) ) then ! Allocate Stuff for advanced LES models 1598 allocate (fwr2(nshg)) 1599 allocate (fwr3(nshg)) 1600 allocate (fwr4(nshg)) 1601 allocate (xavegt(nfath,12)) 1602 allocate (xavegt2(nfath,12)) 1603 allocate (xavegt3(nfath,12)) 1604 allocate (stabdis(nfath)) 1605 endif 1606 1607c.... get dynamic model coefficient 1608c 1609 ilesmod=iLES/10 1610c 1611c digit bit set filter rule, 10 bit set model 1612c 1613 if (ilesmod.eq.0) then ! 0 < iLES< 10 => dyn. model calculated 1614 ! at nodes based on discrete filtering 1615 1616 1617 if(isubmod.eq.2) then 1618 call SUPGdis(y, ac, shgl, shp, 1619 & iper, ilwork, 1620 & nsons, ifath, x, 1621 & iBC, BC, stabdis, xavegt3) 1622 endif 1623 1624 if( ((isubmod.eq.0).or.(isubmod.eq.2)))then ! If no 1625 ! sub-model 1626 ! or SUPG 1627 ! model wanted 1628 1629 if(i2filt.eq.0)then ! If simple filter 1630 1631 if(modlstats .eq. 0) then ! If no model stats wanted 1632 call getdmc (y, shgl, shp, 1633 & iper, ilwork, nsons, 1634 & ifath, x) 1635 else ! else get model stats 1636 call stdfdmc (y, shgl, shp, 1637 & iper, ilwork, nsons, 1638 & ifath, x) 1639 endif ! end of stats if statement 1640 1641 else ! else if twice filtering 1642 1643 call widefdmc(y, shgl, shp, 1644 & iper, ilwork, nsons, 1645 & ifath, x) 1646 1647 1648 endif ! end of simple filter if statement 1649 1650 endif ! end of SUPG or no sub-model if statement 1651 1652 1653 if( (isubmod.eq.1) ) then ! If DFWR sub-model wanted 1654 call cdelBHsq (y, shgl, shp, 1655 & iper, ilwork, nsons, 1656 & ifath, x, cdelsq1) 1657 call FiltRat (y, shgl, shp, 1658 & iper, ilwork, nsons, 1659 & ifath, x, cdelsq1, 1660 & fwr4, fwr3) 1661 1662 1663 if (i2filt.eq.0) then ! If simple filter wanted 1664 call DFWRsfdmc(y, shgl, shp, 1665 & iper, ilwork, nsons, 1666 & ifath, x, fwr2, fwr3) 1667 else ! else if twice filtering wanted 1668 call DFWRwfdmc(y, shgl, shp, 1669 & iper, ilwork, nsons, 1670 & ifath, x, fwr4, fwr4) 1671 endif ! end of simple filter if statement 1672 1673 endif ! end of DFWR sub-model if statement 1674 1675 if( (isubmod.eq.2) )then ! If SUPG sub-model wanted 1676 call dmcSUPG (y, ac, shgl, 1677 & shp, iper, ilwork, 1678 & nsons, ifath, x, 1679 & iBC, BC, rowp, colm, 1680 & xavegt2, stabdis) 1681 endif 1682 1683 if(idis.eq.1)then ! If SUPG/Model dissipation wanted 1684 call ediss (y, ac, shgl, 1685 & shp, iper, ilwork, 1686 & nsons, ifath, x, 1687 & iBC, BC, xavegt) 1688 endif 1689 1690 endif ! end of ilesmod 1691 1692 if (ilesmod .eq. 1) then ! 10 < iLES < 20 => dynamic-mixed 1693 ! at nodes based on discrete filtering 1694 call bardmc (y, shgl, shp, 1695 & iper, ilwork, 1696 & nsons, ifath, x) 1697 endif 1698 1699 if (ilesmod .eq. 2) then ! 20 < iLES < 30 => dynamic at quad 1700 ! pts based on lumped projection filt. 1701 1702 if(isubmod.eq.0)then 1703 call projdmc (y, shgl, shp, 1704 & iper, ilwork, x) 1705 else 1706 call cpjdmcnoi (y, shgl, shp, 1707 & iper, ilwork, x, 1708 & rowp, colm, 1709 & iBC, BC) 1710 endif 1711 1712 endif 1713 1714 if( (iLES.gt.1) ) then ! Deallocate Stuff for advanced LES models 1715 deallocate (fwr2) 1716 deallocate (fwr3) 1717 deallocate (fwr4) 1718 deallocate (xavegt) 1719 deallocate (xavegt2) 1720 deallocate (xavegt3) 1721 deallocate (stabdis) 1722 endif 1723 return 1724 end 1725 1726c 1727c...initialize the coefficients for the impedance convolution 1728c 1729 subroutine CalcImpConvCoef (numISrfs, numTpoints) 1730 1731 use convolImpFlow !uses flow history and impedance for convolution 1732 1733 include "common.h" !for alfi 1734 1735 integer numISrfs, numTpoints 1736 1737 allocate (ConvCoef(numTpoints+2,3)) !same time discret. for all imp. BC 1738 do j=1,numTpoints+2 1739 ConvCoef(j,:)=0.5/numTpoints !dt/2 divided by period T=N*dt 1740 ConvCoef(j,1)=ConvCoef(j,1)*(1.0-alfi)*(1.0-alfi) 1741 ConvCoef(j,2)=ConvCoef(j,2)*(1.0+2*alfi*(1.0-alfi)) 1742 ConvCoef(j,3)=ConvCoef(j,3)*alfi*alfi 1743 enddo 1744 ConvCoef(1,2)=zero 1745 ConvCoef(1,3)=zero 1746 ConvCoef(2,3)=zero 1747 ConvCoef(numTpoints+1,1)=zero 1748 ConvCoef(numTpoints+2,2)=zero 1749 ConvCoef(numTpoints+2,1)=zero 1750c 1751c...calculate the coefficients for the impedance convolution 1752c 1753 allocate (ImpConvCoef(numTpoints+2,numISrfs)) 1754 1755c..coefficients below assume Q linear in time step, Z constant 1756c do j=3,numTpoints 1757c ImpConvCoef(j,:) = ValueListImp(j-1,:)*ConvCoef(j,3) 1758c & + ValueListImp(j,:)*ConvCoef(j,2) 1759c & + ValueListImp(j+1,:)*ConvCoef(j,1) 1760c enddo 1761c ImpConvCoef(1,:) = ValueListImp(2,:)*ConvCoef(1,1) 1762c ImpConvCoef(2,:) = ValueListImp(2,:)*ConvCoef(2,2) 1763c & + ValueListImp(3,:)*ConvCoef(2,1) 1764c ImpConvCoef(numTpoints+1,:) = 1765c & ValueListImp(numTpoints,:)*ConvCoef(numTpoints+1,3) 1766c & + ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+1,2) 1767c ImpConvCoef(numTpoints+2,:) = 1768c & ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+2,3) 1769 1770c..try easiest convolution Q and Z constant per time step 1771 do j=3,numTpoints+1 1772 ImpConvCoef(j,:) = ValueListImp(j-1,:)/numTpoints 1773 enddo 1774 ImpConvCoef(1,:) =zero 1775 ImpConvCoef(2,:) =zero 1776 ImpConvCoef(numTpoints+2,:) = 1777 & ValueListImp(numTpoints+1,:)/numTpoints 1778c compensate for yalpha passed not y in Elmgmr() 1779 ImpConvCoef(numTpoints+1,:)= ImpConvCoef(numTpoints+1,:) 1780 & - ImpConvCoef(numTpoints+2,:)*(1.0-alfi)/alfi 1781 ImpConvCoef(numTpoints+2,:)= ImpConvCoef(numTpoints+2,:)/alfi 1782 return 1783 end 1784 1785c 1786c ... update the flow rate history for the impedance convolution, filter it and write it out 1787c 1788 subroutine UpdHistConv(y,nsrfIdList,numSrfs) 1789 1790 use convolImpFlow !brings ntimeptpT, QHistImp, QHistTry, QHistTryF, numImpSrfs 1791 use convolRCRFlow !brings QHistRCR, numRCRSrfs 1792 1793 include "common.h" 1794 1795 integer nsrfIdList(0:MAXSURF), numSrfs 1796 character*20 fname1 1797 character*10 cname2 1798 character*5 cname 1799 real*8 y(nshg,3) !velocity at time n+1 1800 real*8 NewQ(0:MAXSURF) !temporary unknown for the flow rate 1801 !that needs to be added to the flow history 1802 1803 call GetFlowQ(NewQ,y,nsrfIdList,numSrfs) !new flow at time n+1 1804c 1805c... for imp BC: shift QHist, add new constribution, filter and write out 1806c 1807 if(numImpSrfs.gt.zero) then 1808 do j=1, ntimeptpT 1809 QHistImp(j,1:numSrfs)=QHistImp(j+1,1:numSrfs) 1810 enddo 1811 QHistImp(ntimeptpT+1,1:numSrfs) = NewQ(1:numSrfs) 1812 1813c 1814c....filter the flow rate history 1815c 1816 cutfreq = 10 !hardcoded cutting frequency of the filter 1817 do j=1, ntimeptpT 1818 QHistTry(j,:)=QHistImp(j+1,:) 1819 enddo 1820 call Filter(QHistTryF,QHistTry,ntimeptpT,Delt(1),cutfreq) 1821c.... if no filter applied then uncomment next three lines 1822c do j=1, ntimeptpT 1823c QHistTryF(j,:)=QHistTry(j,:) 1824c enddo 1825 1826c QHistImp(1,:)=zero ! why do we do this? for beta(1,:) = zero it does not really matters 1827 do j=1, ntimeptpT 1828 QHistImp(j+1,:)=QHistTryF(j,:) 1829 enddo 1830c 1831c.... write out the new history of flow rates to Qhistor.dat 1832c 1833 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1834 & (istep .eq. nstep(1)))) .and. 1835 & (myrank .eq. master)) then 1836 open(unit=816, file='Qhistor.dat',status='replace') 1837 write(816,*) ntimeptpT 1838 do j=1,ntimeptpT+1 1839 write(816,*) (QHistImp(j,n),n=1, numSrfs) 1840 enddo 1841 close(816) 1842c... write out a copy with step number to be able to restart 1843 fname1 = 'Qhistor' 1844 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1845 open(unit=8166,file=trim(fname1),status='unknown') 1846 write(8166,*) ntimeptpT 1847 do j=1,ntimeptpT+1 1848 write(8166,*) (QHistImp(j,n),n=1, numSrfs) 1849 enddo 1850 close(8166) 1851 endif 1852 endif 1853 1854c 1855c... for RCR bc just add the new contribution 1856c 1857 if(numRCRSrfs.gt.zero) then 1858 QHistRCR(lstep+1,1:numSrfs) = NewQ(1:numSrfs) 1859c 1860c.... write out the new history of flow rates to Qhistor.dat 1861c 1862 if ((irs .ge. 1) .and. (myrank .eq. master)) then 1863 if(istep.eq.1) then 1864 open(unit=816,file='Qhistor.dat',status='unknown') 1865 else 1866 open(unit=816,file='Qhistor.dat',position='append') 1867 endif 1868 if(istep.eq.1) then 1869 do j=1,lstep 1870 write(816,*) j, (QHistRCR(j,n),n=1,numSrfs) ! read from file of previous run 1871 enddo 1872 endif 1873 write(816,*) lstep+1, (QHistRCR(lstep+1,n),n=1, numSrfs) 1874 close(816) 1875c... write out a copy with step number to be able to restart 1876 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1877 & (istep .eq. nstep(1)))) .and. 1878 & (myrank .eq. master)) then 1879 fname1 = 'Qhistor' 1880 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1881 open(unit=8166,file=trim(fname1),status='unknown') 1882 write(8166,*) lstep+1 1883 do j=1,lstep+1 1884 write(8166,*) (QHistRCR(j,n),n=1, numSrfs) 1885 enddo 1886 close(8166) 1887 endif 1888 endif 1889 endif 1890 1891 return 1892 end 1893 1894c 1895c...calculate the time varying coefficients for the RCR convolution 1896c 1897 subroutine CalcRCRConvCoef (stepn, numSrfs) 1898 1899 use convolRCRFlow !brings in ValueListRCR, dtRCR 1900 1901 include "common.h" !brings alfi 1902 1903 integer numSrfs, stepn 1904 1905 RCRConvCoef = zero 1906 if (stepn .eq. 0) then 1907 RCRConvCoef(1,:) = ValueListRCR(1,:)*(1.0-alfi) + 1908 & ValueListRCR(3,:)*(-alfi + 1.0 + 1/dtRCR(:) 1909 & - exp(-alfi*dtRCR(:))*(1 + 1/dtRCR(:))) 1910 RCRConvCoef(2,:) = ValueListRCR(1,:)*alfi 1911 & + ValueListRCR(3,:) 1912 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1913 endif 1914 if (stepn .ge. 1) then 1915 RCRConvCoef(1,:) =-ValueListRCR(3,:)*exp(-dtRCR(:)*(stepn+alfi)) 1916 & *(1 + (1 - exp(dtRCR(:)))/dtRCR(:)) 1917 RCRConvCoef(stepn+1,:) = ValueListRCR(1,:)*(1-alfi) 1918 & - ValueListRCR(3,:)*(alfi - 1 - 1/dtRCR(:) 1919 & + exp(-alfi*dtRCR(:))/dtRCR(:)*(2 - exp(-dtRCR(:)))) 1920 RCRConvCoef(stepn+2,:) = ValueListRCR(1,:)*alfi 1921 & + ValueListRCR(3,:) 1922 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1923 endif 1924 if (stepn .ge. 2) then 1925 do j=2,stepn 1926 RCRConvCoef(j,:) = ValueListRCR(3,:)/dtRCR(:)* 1927 & exp(-dtRCR(:)*(stepn + alfi + 2 - j))* 1928 & (1 - exp(dtRCR(:)))**2 1929 enddo 1930 endif 1931 1932c compensate for yalpha passed not y in Elmgmr() 1933 RCRConvCoef(stepn+1,:)= RCRConvCoef(stepn+1,:) 1934 & - RCRConvCoef(stepn+2,:)*(1.0-alfi)/alfi 1935 RCRConvCoef(stepn+2,:)= RCRConvCoef(stepn+2,:)/alfi 1936 1937 return 1938 end 1939 1940c 1941c...calculate the time dependent H operator for the RCR convolution 1942c 1943 subroutine CalcHopRCR (timestepRCR, stepn, numSrfs) 1944 1945 use convolRCRFlow !brings in HopRCR, dtRCR 1946 1947 include "common.h" 1948 1949 integer numSrfs, stepn 1950 real*8 PdistCur(0:MAXSURF), timestepRCR 1951 1952 HopRCR=zero 1953 call RCRint(timestepRCR*(stepn + alfi),PdistCur) 1954 HopRCR(1:numSrfs) = RCRic(1:numSrfs) 1955 & *exp(-dtRCR(1:numSrfs)*(stepn + alfi)) + PdistCur(1:numSrfs) 1956 return 1957 end 1958c 1959c ... initialize the influence of the initial conditions for the RCR BC 1960c 1961 subroutine calcRCRic(y,srfIdList,numSrfs) 1962 1963 use convolRCRFlow !brings RCRic, ValueListRCR, ValuePdist 1964 1965 include "common.h" 1966 1967 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled 1968 real*8 y(nshg,4) !need velocity and pressure 1969 real*8 Qini(0:MAXSURF) !initial flow rate 1970 real*8 PdistIni(0:MAXSURF) !initial distal pressure 1971 real*8 Pini(0:MAXSURF),CoupleArea(0:MAXSURF) ! initial pressure 1972 real*8 VelOnly(nshg,3), POnly(nshg) 1973 1974 allocate (RCRic(0:MAXSURF)) 1975 1976 if(lstep.eq.0) then 1977 VelOnly(:,1:3)=y(:,1:3) 1978 call GetFlowQ(Qini,VelOnly,srfIdList,numSrfs) !get initial flow 1979 QHistRCR(1,1:numSrfs)=Qini(1:numSrfs) !initialize QHistRCR 1980 POnly(:)=y(:,4) ! pressure 1981 call integrScalar(Pini,POnly,srfIdList,numSrfs) !get initial pressure integral 1982 POnly(:)=one ! one to get area 1983 call integrScalar(CoupleArea,POnly,srfIdList,numSrfs) !get surf area 1984 Pini(1:numSrfs) = Pini(1:numSrfs)/CoupleArea(1:numSrfs) 1985 else 1986 Qini(1:numSrfs)=QHistRCR(1,1:numSrfs) 1987 Pini(1:numSrfs)=zero ! hack 1988 endif 1989 call RCRint(istep,PdistIni) !get initial distal P (use istep) 1990 RCRic(1:numSrfs) = Pini(1:numSrfs) 1991 & - ValueListRCR(1,:)*Qini(1:numSrfs)-PdistIni(1:numSrfs) 1992 return 1993 end 1994 1995c.........function that integrates a scalar over a boundary 1996 subroutine integrScalar(scalInt,scal,srfIdList,numSrfs) 1997 1998 use pvsQbi !brings ndsurf, NASC 1999 2000 include "common.h" 2001 include "mpif.h" 2002 2003 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled, i, k 2004 real*8 scal(nshg), scalInt(0:MAXSURF), scalIntProc(0:MAXSURF) 2005 2006 scalIntProc = zero 2007 do i = 1,nshg 2008 if(numSrfs.gt.zero) then 2009 do k = 1,numSrfs 2010 irankCoupled = 0 2011 if (srfIdList(k).eq.ndsurf(i)) then 2012 irankCoupled=k 2013 scalIntProc(irankCoupled) = scalIntProc(irankCoupled) 2014 & + NASC(i)*scal(i) 2015 exit 2016 endif 2017 enddo 2018 endif 2019 enddo 2020c 2021c at this point, each scalint has its "nodes" contributions to the scalar 2022c accumulated into scalIntProc. Note, because NASC is on processor this 2023c will NOT be the scalar for the surface yet 2024c 2025c.... reduce integrated scalar for each surface, push on scalInt 2026c 2027 npars=MAXSURF+1 2028 call MPI_ALLREDUCE (scalIntProc, scalInt(:), npars, 2029 & MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr) 2030 2031 return 2032 end 2033 2034 subroutine writeTimingMessage(key,iomode,timing) 2035 use iso_c_binding 2036 use phstr 2037 implicit none 2038 2039 character(len=*) :: key 2040 integer :: iomode 2041 real*8 :: timing 2042 character(len=1024) :: timing_msg 2043 character(len=*), parameter :: 2044 & streamModeString = c_char_"stream"//c_null_char, 2045 & fileModeString = c_char_"disk"//c_null_char 2046 2047 timing_msg = c_char_"Time to write "//c_null_char 2048 call phstr_appendStr(timing_msg,key) 2049 if ( iomode .eq. -1 ) then 2050 call phstr_appendStr(timing_msg, streamModeString) 2051 else 2052 call phstr_appendStr(timing_msg, fileModeString) 2053 endif 2054 call phstr_appendStr(timing_msg, c_char_' = '//c_null_char) 2055 call phstr_appendDbl(timing_msg, timing) 2056 write(6,*) trim(timing_msg) 2057 return 2058 end subroutine 2059 2060