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