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