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 1485 3000 continue 1486 1487 1488c 1489c.... close history and aerodynamic forces files 1490c 1491 if (myrank .eq. master) then 1492! close (ihist) 1493 close (iforce) 1494 close(76) 1495 if(numImpSrfs.gt.0 .or. numRCRSrfs.gt.0) then 1496 close (8177) 1497 endif 1498 endif 1499c 1500c.... close varts file for probes 1501c 1502 if(exts) then 1503 do jj=1,ntspts 1504 if (myrank == iv_rank(jj)) then 1505 close(1000+jj) 1506 endif 1507 enddo 1508 deallocate (ivarts) 1509 deallocate (ivartsg) 1510 deallocate (iv_rank) 1511 deallocate (vartssoln) 1512 deallocate (vartssolng) 1513 endif 1514 1515!MR CHANGE 1516 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1517 if(myrank.eq.0) then 1518 write(*,*) 'itrdrv - done with aerodynamic forces' 1519 endif 1520!MR CHANGE 1521 1522 do isrf = 0,MAXSURF 1523! if ( nsrflist(isrf).ne.0 ) then 1524 if ( nsrflist(isrf).ne.0 .and. 1525 & myrank.eq.irankfilesforce(isrf)) then 1526 iunit=60+isrf 1527 close(iunit) 1528 endif 1529 enddo 1530 1531!MR CHANGE 1532 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1533 if(myrank.eq.0) then 1534 write(*,*) 'itrdrv - done with MAXSURF' 1535 endif 1536!MR CHANGE 1537 1538 1539 5 format(1X,F15.10,3X,F15.10,3X,F15.10,3X,F15.10) 1540 444 format(6(2x,e14.7)) 1541c 1542c.... end 1543c 1544 if(nsolflow.eq.1) then 1545 deallocate (lhsK) 1546 deallocate (lhsP) 1547 IF (svLSFlag .NE. 1) THEN 1548 deallocate (aperm) 1549 deallocate (atemp) 1550 ENDIF 1551 endif 1552 if(nsclrsol.gt.0) then 1553 deallocate (lhsS) 1554 deallocate (apermS) 1555 deallocate (atempS) 1556 endif 1557 1558 if(iabc==1) deallocate(acs) 1559 1560!MR CHANGE 1561 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1562 if(myrank.eq.0) then 1563 write(*,*) 'itrdrv - done - BACK TO process.f' 1564 endif 1565!MR CHANGE 1566 1567 1568 1569 return 1570 end 1571 1572 subroutine lesmodels(y, ac, shgl, shp, 1573 & iper, ilwork, rowp, colm, 1574 & nsons, ifath, x, 1575 & iBC, BC) 1576 1577 include "common.h" 1578 1579 real*8 y(nshg,ndof), ac(nshg,ndof), 1580 & x(numnp,nsd), 1581 & BC(nshg,ndofBC) 1582 real*8 shp(MAXTOP,maxsh,MAXQPT), 1583 & shgl(MAXTOP,nsd,maxsh,MAXQPT) 1584 1585c 1586 integer rowp(nshg,nnz), colm(nshg+1), 1587 & iBC(nshg), 1588 & ilwork(nlwork), 1589 & iper(nshg) 1590 dimension ifath(numnp), nsons(nfath) 1591 1592 real*8, allocatable, dimension(:) :: fwr2,fwr3,fwr4 1593 real*8, allocatable, dimension(:) :: stabdis,cdelsq1 1594 real*8, allocatable, dimension(:,:) :: xavegt, xavegt2,xavegt3 1595 1596 if( (iLES.gt.1) ) then ! Allocate Stuff for advanced LES models 1597 allocate (fwr2(nshg)) 1598 allocate (fwr3(nshg)) 1599 allocate (fwr4(nshg)) 1600 allocate (xavegt(nfath,12)) 1601 allocate (xavegt2(nfath,12)) 1602 allocate (xavegt3(nfath,12)) 1603 allocate (stabdis(nfath)) 1604 endif 1605 1606c.... get dynamic model coefficient 1607c 1608 ilesmod=iLES/10 1609c 1610c digit bit set filter rule, 10 bit set model 1611c 1612 if (ilesmod.eq.0) then ! 0 < iLES< 10 => dyn. model calculated 1613 ! at nodes based on discrete filtering 1614 1615 1616 if(isubmod.eq.2) then 1617 call SUPGdis(y, ac, shgl, shp, 1618 & iper, ilwork, 1619 & nsons, ifath, x, 1620 & iBC, BC, stabdis, xavegt3) 1621 endif 1622 1623 if( ((isubmod.eq.0).or.(isubmod.eq.2)))then ! If no 1624 ! sub-model 1625 ! or SUPG 1626 ! model wanted 1627 1628 if(i2filt.eq.0)then ! If simple filter 1629 1630 if(modlstats .eq. 0) then ! If no model stats wanted 1631 call getdmc (y, shgl, shp, 1632 & iper, ilwork, nsons, 1633 & ifath, x) 1634 else ! else get model stats 1635 call stdfdmc (y, shgl, shp, 1636 & iper, ilwork, nsons, 1637 & ifath, x) 1638 endif ! end of stats if statement 1639 1640 else ! else if twice filtering 1641 1642 call widefdmc(y, shgl, shp, 1643 & iper, ilwork, nsons, 1644 & ifath, x) 1645 1646 1647 endif ! end of simple filter if statement 1648 1649 endif ! end of SUPG or no sub-model if statement 1650 1651 1652 if( (isubmod.eq.1) ) then ! If DFWR sub-model wanted 1653 call cdelBHsq (y, shgl, shp, 1654 & iper, ilwork, nsons, 1655 & ifath, x, cdelsq1) 1656 call FiltRat (y, shgl, shp, 1657 & iper, ilwork, nsons, 1658 & ifath, x, cdelsq1, 1659 & fwr4, fwr3) 1660 1661 1662 if (i2filt.eq.0) then ! If simple filter wanted 1663 call DFWRsfdmc(y, shgl, shp, 1664 & iper, ilwork, nsons, 1665 & ifath, x, fwr2, fwr3) 1666 else ! else if twice filtering wanted 1667 call DFWRwfdmc(y, shgl, shp, 1668 & iper, ilwork, nsons, 1669 & ifath, x, fwr4, fwr4) 1670 endif ! end of simple filter if statement 1671 1672 endif ! end of DFWR sub-model if statement 1673 1674 if( (isubmod.eq.2) )then ! If SUPG sub-model wanted 1675 call dmcSUPG (y, ac, shgl, 1676 & shp, iper, ilwork, 1677 & nsons, ifath, x, 1678 & iBC, BC, rowp, colm, 1679 & xavegt2, stabdis) 1680 endif 1681 1682 if(idis.eq.1)then ! If SUPG/Model dissipation wanted 1683 call ediss (y, ac, shgl, 1684 & shp, iper, ilwork, 1685 & nsons, ifath, x, 1686 & iBC, BC, xavegt) 1687 endif 1688 1689 endif ! end of ilesmod 1690 1691 if (ilesmod .eq. 1) then ! 10 < iLES < 20 => dynamic-mixed 1692 ! at nodes based on discrete filtering 1693 call bardmc (y, shgl, shp, 1694 & iper, ilwork, 1695 & nsons, ifath, x) 1696 endif 1697 1698 if (ilesmod .eq. 2) then ! 20 < iLES < 30 => dynamic at quad 1699 ! pts based on lumped projection filt. 1700 1701 if(isubmod.eq.0)then 1702 call projdmc (y, shgl, shp, 1703 & iper, ilwork, x) 1704 else 1705 call cpjdmcnoi (y, shgl, shp, 1706 & iper, ilwork, x, 1707 & rowp, colm, 1708 & iBC, BC) 1709 endif 1710 1711 endif 1712 1713 if( (iLES.gt.1) ) then ! Deallocate Stuff for advanced LES models 1714 deallocate (fwr2) 1715 deallocate (fwr3) 1716 deallocate (fwr4) 1717 deallocate (xavegt) 1718 deallocate (xavegt2) 1719 deallocate (xavegt3) 1720 deallocate (stabdis) 1721 endif 1722 return 1723 end 1724 1725c 1726c...initialize the coefficients for the impedance convolution 1727c 1728 subroutine CalcImpConvCoef (numISrfs, numTpoints) 1729 1730 use convolImpFlow !uses flow history and impedance for convolution 1731 1732 include "common.h" !for alfi 1733 1734 integer numISrfs, numTpoints 1735 1736 allocate (ConvCoef(numTpoints+2,3)) !same time discret. for all imp. BC 1737 do j=1,numTpoints+2 1738 ConvCoef(j,:)=0.5/numTpoints !dt/2 divided by period T=N*dt 1739 ConvCoef(j,1)=ConvCoef(j,1)*(1.0-alfi)*(1.0-alfi) 1740 ConvCoef(j,2)=ConvCoef(j,2)*(1.0+2*alfi*(1.0-alfi)) 1741 ConvCoef(j,3)=ConvCoef(j,3)*alfi*alfi 1742 enddo 1743 ConvCoef(1,2)=zero 1744 ConvCoef(1,3)=zero 1745 ConvCoef(2,3)=zero 1746 ConvCoef(numTpoints+1,1)=zero 1747 ConvCoef(numTpoints+2,2)=zero 1748 ConvCoef(numTpoints+2,1)=zero 1749c 1750c...calculate the coefficients for the impedance convolution 1751c 1752 allocate (ImpConvCoef(numTpoints+2,numISrfs)) 1753 1754c..coefficients below assume Q linear in time step, Z constant 1755c do j=3,numTpoints 1756c ImpConvCoef(j,:) = ValueListImp(j-1,:)*ConvCoef(j,3) 1757c & + ValueListImp(j,:)*ConvCoef(j,2) 1758c & + ValueListImp(j+1,:)*ConvCoef(j,1) 1759c enddo 1760c ImpConvCoef(1,:) = ValueListImp(2,:)*ConvCoef(1,1) 1761c ImpConvCoef(2,:) = ValueListImp(2,:)*ConvCoef(2,2) 1762c & + ValueListImp(3,:)*ConvCoef(2,1) 1763c ImpConvCoef(numTpoints+1,:) = 1764c & ValueListImp(numTpoints,:)*ConvCoef(numTpoints+1,3) 1765c & + ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+1,2) 1766c ImpConvCoef(numTpoints+2,:) = 1767c & ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+2,3) 1768 1769c..try easiest convolution Q and Z constant per time step 1770 do j=3,numTpoints+1 1771 ImpConvCoef(j,:) = ValueListImp(j-1,:)/numTpoints 1772 enddo 1773 ImpConvCoef(1,:) =zero 1774 ImpConvCoef(2,:) =zero 1775 ImpConvCoef(numTpoints+2,:) = 1776 & ValueListImp(numTpoints+1,:)/numTpoints 1777c compensate for yalpha passed not y in Elmgmr() 1778 ImpConvCoef(numTpoints+1,:)= ImpConvCoef(numTpoints+1,:) 1779 & - ImpConvCoef(numTpoints+2,:)*(1.0-alfi)/alfi 1780 ImpConvCoef(numTpoints+2,:)= ImpConvCoef(numTpoints+2,:)/alfi 1781 return 1782 end 1783 1784c 1785c ... update the flow rate history for the impedance convolution, filter it and write it out 1786c 1787 subroutine UpdHistConv(y,nsrfIdList,numSrfs) 1788 1789 use convolImpFlow !brings ntimeptpT, QHistImp, QHistTry, QHistTryF, numImpSrfs 1790 use convolRCRFlow !brings QHistRCR, numRCRSrfs 1791 1792 include "common.h" 1793 1794 integer nsrfIdList(0:MAXSURF), numSrfs 1795 character*20 fname1 1796 character*10 cname2 1797 character*5 cname 1798 real*8 y(nshg,3) !velocity at time n+1 1799 real*8 NewQ(0:MAXSURF) !temporary unknown for the flow rate 1800 !that needs to be added to the flow history 1801 1802 call GetFlowQ(NewQ,y,nsrfIdList,numSrfs) !new flow at time n+1 1803c 1804c... for imp BC: shift QHist, add new constribution, filter and write out 1805c 1806 if(numImpSrfs.gt.zero) then 1807 do j=1, ntimeptpT 1808 QHistImp(j,1:numSrfs)=QHistImp(j+1,1:numSrfs) 1809 enddo 1810 QHistImp(ntimeptpT+1,1:numSrfs) = NewQ(1:numSrfs) 1811 1812c 1813c....filter the flow rate history 1814c 1815 cutfreq = 10 !hardcoded cutting frequency of the filter 1816 do j=1, ntimeptpT 1817 QHistTry(j,:)=QHistImp(j+1,:) 1818 enddo 1819 call Filter(QHistTryF,QHistTry,ntimeptpT,Delt(1),cutfreq) 1820c.... if no filter applied then uncomment next three lines 1821c do j=1, ntimeptpT 1822c QHistTryF(j,:)=QHistTry(j,:) 1823c enddo 1824 1825c QHistImp(1,:)=zero ! why do we do this? for beta(1,:) = zero it does not really matters 1826 do j=1, ntimeptpT 1827 QHistImp(j+1,:)=QHistTryF(j,:) 1828 enddo 1829c 1830c.... write out the new history of flow rates to Qhistor.dat 1831c 1832 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1833 & (istep .eq. nstep(1)))) .and. 1834 & (myrank .eq. master)) then 1835 open(unit=816, file='Qhistor.dat',status='replace') 1836 write(816,*) ntimeptpT 1837 do j=1,ntimeptpT+1 1838 write(816,*) (QHistImp(j,n),n=1, numSrfs) 1839 enddo 1840 close(816) 1841c... write out a copy with step number to be able to restart 1842 fname1 = 'Qhistor' 1843 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1844 open(unit=8166,file=trim(fname1),status='unknown') 1845 write(8166,*) ntimeptpT 1846 do j=1,ntimeptpT+1 1847 write(8166,*) (QHistImp(j,n),n=1, numSrfs) 1848 enddo 1849 close(8166) 1850 endif 1851 endif 1852 1853c 1854c... for RCR bc just add the new contribution 1855c 1856 if(numRCRSrfs.gt.zero) then 1857 QHistRCR(lstep+1,1:numSrfs) = NewQ(1:numSrfs) 1858c 1859c.... write out the new history of flow rates to Qhistor.dat 1860c 1861 if ((irs .ge. 1) .and. (myrank .eq. master)) then 1862 if(istep.eq.1) then 1863 open(unit=816,file='Qhistor.dat',status='unknown') 1864 else 1865 open(unit=816,file='Qhistor.dat',position='append') 1866 endif 1867 if(istep.eq.1) then 1868 do j=1,lstep 1869 write(816,*) j, (QHistRCR(j,n),n=1,numSrfs) ! read from file of previous run 1870 enddo 1871 endif 1872 write(816,*) lstep+1, (QHistRCR(lstep+1,n),n=1, numSrfs) 1873 close(816) 1874c... write out a copy with step number to be able to restart 1875 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1876 & (istep .eq. nstep(1)))) .and. 1877 & (myrank .eq. master)) then 1878 fname1 = 'Qhistor' 1879 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1880 open(unit=8166,file=trim(fname1),status='unknown') 1881 write(8166,*) lstep+1 1882 do j=1,lstep+1 1883 write(8166,*) (QHistRCR(j,n),n=1, numSrfs) 1884 enddo 1885 close(8166) 1886 endif 1887 endif 1888 endif 1889 1890 return 1891 end 1892 1893c 1894c...calculate the time varying coefficients for the RCR convolution 1895c 1896 subroutine CalcRCRConvCoef (stepn, numSrfs) 1897 1898 use convolRCRFlow !brings in ValueListRCR, dtRCR 1899 1900 include "common.h" !brings alfi 1901 1902 integer numSrfs, stepn 1903 1904 RCRConvCoef = zero 1905 if (stepn .eq. 0) then 1906 RCRConvCoef(1,:) = ValueListRCR(1,:)*(1.0-alfi) + 1907 & ValueListRCR(3,:)*(-alfi + 1.0 + 1/dtRCR(:) 1908 & - exp(-alfi*dtRCR(:))*(1 + 1/dtRCR(:))) 1909 RCRConvCoef(2,:) = ValueListRCR(1,:)*alfi 1910 & + ValueListRCR(3,:) 1911 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1912 endif 1913 if (stepn .ge. 1) then 1914 RCRConvCoef(1,:) =-ValueListRCR(3,:)*exp(-dtRCR(:)*(stepn+alfi)) 1915 & *(1 + (1 - exp(dtRCR(:)))/dtRCR(:)) 1916 RCRConvCoef(stepn+1,:) = ValueListRCR(1,:)*(1-alfi) 1917 & - ValueListRCR(3,:)*(alfi - 1 - 1/dtRCR(:) 1918 & + exp(-alfi*dtRCR(:))/dtRCR(:)*(2 - exp(-dtRCR(:)))) 1919 RCRConvCoef(stepn+2,:) = ValueListRCR(1,:)*alfi 1920 & + ValueListRCR(3,:) 1921 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1922 endif 1923 if (stepn .ge. 2) then 1924 do j=2,stepn 1925 RCRConvCoef(j,:) = ValueListRCR(3,:)/dtRCR(:)* 1926 & exp(-dtRCR(:)*(stepn + alfi + 2 - j))* 1927 & (1 - exp(dtRCR(:)))**2 1928 enddo 1929 endif 1930 1931c compensate for yalpha passed not y in Elmgmr() 1932 RCRConvCoef(stepn+1,:)= RCRConvCoef(stepn+1,:) 1933 & - RCRConvCoef(stepn+2,:)*(1.0-alfi)/alfi 1934 RCRConvCoef(stepn+2,:)= RCRConvCoef(stepn+2,:)/alfi 1935 1936 return 1937 end 1938 1939c 1940c...calculate the time dependent H operator for the RCR convolution 1941c 1942 subroutine CalcHopRCR (timestepRCR, stepn, numSrfs) 1943 1944 use convolRCRFlow !brings in HopRCR, dtRCR 1945 1946 include "common.h" 1947 1948 integer numSrfs, stepn 1949 real*8 PdistCur(0:MAXSURF), timestepRCR 1950 1951 HopRCR=zero 1952 call RCRint(timestepRCR*(stepn + alfi),PdistCur) 1953 HopRCR(1:numSrfs) = RCRic(1:numSrfs) 1954 & *exp(-dtRCR(1:numSrfs)*(stepn + alfi)) + PdistCur(1:numSrfs) 1955 return 1956 end 1957c 1958c ... initialize the influence of the initial conditions for the RCR BC 1959c 1960 subroutine calcRCRic(y,srfIdList,numSrfs) 1961 1962 use convolRCRFlow !brings RCRic, ValueListRCR, ValuePdist 1963 1964 include "common.h" 1965 1966 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled 1967 real*8 y(nshg,4) !need velocity and pressure 1968 real*8 Qini(0:MAXSURF) !initial flow rate 1969 real*8 PdistIni(0:MAXSURF) !initial distal pressure 1970 real*8 Pini(0:MAXSURF),CoupleArea(0:MAXSURF) ! initial pressure 1971 real*8 VelOnly(nshg,3), POnly(nshg) 1972 1973 allocate (RCRic(0:MAXSURF)) 1974 1975 if(lstep.eq.0) then 1976 VelOnly(:,1:3)=y(:,1:3) 1977 call GetFlowQ(Qini,VelOnly,srfIdList,numSrfs) !get initial flow 1978 QHistRCR(1,1:numSrfs)=Qini(1:numSrfs) !initialize QHistRCR 1979 POnly(:)=y(:,4) ! pressure 1980 call integrScalar(Pini,POnly,srfIdList,numSrfs) !get initial pressure integral 1981 POnly(:)=one ! one to get area 1982 call integrScalar(CoupleArea,POnly,srfIdList,numSrfs) !get surf area 1983 Pini(1:numSrfs) = Pini(1:numSrfs)/CoupleArea(1:numSrfs) 1984 else 1985 Qini(1:numSrfs)=QHistRCR(1,1:numSrfs) 1986 Pini(1:numSrfs)=zero ! hack 1987 endif 1988 call RCRint(istep,PdistIni) !get initial distal P (use istep) 1989 RCRic(1:numSrfs) = Pini(1:numSrfs) 1990 & - ValueListRCR(1,:)*Qini(1:numSrfs)-PdistIni(1:numSrfs) 1991 return 1992 end 1993 1994c.........function that integrates a scalar over a boundary 1995 subroutine integrScalar(scalInt,scal,srfIdList,numSrfs) 1996 1997 use pvsQbi !brings ndsurf, NASC 1998 1999 include "common.h" 2000 include "mpif.h" 2001 2002 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled, i, k 2003 real*8 scal(nshg), scalInt(0:MAXSURF), scalIntProc(0:MAXSURF) 2004 2005 scalIntProc = zero 2006 do i = 1,nshg 2007 if(numSrfs.gt.zero) then 2008 do k = 1,numSrfs 2009 irankCoupled = 0 2010 if (srfIdList(k).eq.ndsurf(i)) then 2011 irankCoupled=k 2012 scalIntProc(irankCoupled) = scalIntProc(irankCoupled) 2013 & + NASC(i)*scal(i) 2014 exit 2015 endif 2016 enddo 2017 endif 2018 enddo 2019c 2020c at this point, each scalint has its "nodes" contributions to the scalar 2021c accumulated into scalIntProc. Note, because NASC is on processor this 2022c will NOT be the scalar for the surface yet 2023c 2024c.... reduce integrated scalar for each surface, push on scalInt 2025c 2026 npars=MAXSURF+1 2027 call MPI_ALLREDUCE (scalIntProc, scalInt(:), npars, 2028 & MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr) 2029 2030 return 2031 end 2032 2033 subroutine writeTimingMessage(key,iomode,timing) 2034 use iso_c_binding 2035 use phstr 2036 implicit none 2037 2038 character(len=*) :: key 2039 integer :: iomode 2040 real*8 :: timing 2041 character(len=1024) :: timing_msg 2042 character(len=*), parameter :: 2043 & streamModeString = c_char_"stream"//c_null_char, 2044 & fileModeString = c_char_"disk"//c_null_char 2045 2046 timing_msg = c_char_"Time to write "//c_null_char 2047 call phstr_appendStr(timing_msg,key) 2048 if ( iomode .eq. -1 ) then 2049 call phstr_appendStr(timing_msg, streamModeString) 2050 else 2051 call phstr_appendStr(timing_msg, fileModeString) 2052 endif 2053 call phstr_appendStr(timing_msg, c_char_' = '//c_null_char) 2054 call phstr_appendDbl(timing_msg, timing) 2055 write(6,*) trim(timing_msg) 2056 return 2057 end subroutine 2058 2059