subroutine itrdrv (y, ac, uold, x, & iBC, BC, & iper, ilwork, shp, & shgl, shpb, shglb, & ifath, velbar, nsons ) c c---------------------------------------------------------------------- c c This iterative driver is the semi-discrete, predictor multi-corrector c algorithm. It contains the Hulbert Generalized Alpha method which c is 2nd order accurate for Rho_inf from 0 to 1. The method can be c made first-order accurate by setting Rho_inf=-1. It uses a c GMRES iterative solver. c c working arrays: c y (nshg,ndof) : Y variables c x (nshg,nsd) : node coordinates c iBC (nshg) : BC codes c BC (nshg,ndofBC) : BC constraint parameters c iper (nshg) : periodicity table c c shape functions: c shp (nshape,ngauss) : interior element shape functions c shgl (nsd,nshape,ngauss) : local shape function gradients c shpb (nshapeb,ngaussb) : boundary element shape functions c shglb (nsd,nshapeb,ngaussb) : bdry. elt. shape gradients c c Zdenek Johan, Winter 1991. (Fortran 90) c---------------------------------------------------------------------- c use pvsQbi !gives us splag (the spmass at the end of this run use specialBC !gives us itvn use timedata !allows collection of time series use turbSA include "common.h" include "mpif.h" include "auxmpi.h" c dimension y(nshg,ndof), ac(nshg,ndof), & yold(nshg,ndof), acold(nshg,ndof), & x(numnp,nsd), iBC(nshg), & BC(nshg,ndofBC), ilwork(nlwork), & iper(nshg), uold(nshg,nsd) c dimension res(nshg,nflow), BDiag(nshg,nflow,nflow), & rest(nshg), solinc(nshg,ndof) c dimension shp(MAXTOP,maxsh,MAXQPT), & shgl(MAXTOP,nsd,maxsh,MAXQPT), & shpb(MAXTOP,maxsh,MAXQPT), & shglb(MAXTOP,nsd,maxsh,MAXQPT) real*8 almit, alfit, gamit dimension ifath(numnp), velbar(nfath,ndof), nsons(nfath) real*8 rerr(nshg,10),ybar(nshg,ndof+8) ! 8 is for avg. of square as uu, vv, ww, pp, TT, uv, uw, and vw integer, allocatable, dimension(:) :: ivarts integer, allocatable, dimension(:) :: ivartsg real*8, allocatable, dimension(:) :: vartssoln real*8, allocatable, dimension(:) :: vartssolng real*8, allocatable, dimension(:,:,:) :: vartsbuff ! assuming three profiles to control (inlet, bottom FC and top FC) ! store velocity profile set via BC real*8 vbc_prof(nshg,3) character*20 fname1, fmt1, fname2, fmt2 character*4 fname4c ! 4 characters character*60 fvarts character*5 cname character*10 cname2 integer ifuncs(6), iarray(10) real*8 elDw(numel) ! element average of DES d variable c c Here are the data structures for sparse matrix GMRES c integer, allocatable, dimension(:,:) :: rowp integer, allocatable, dimension(:) :: colm real*8, allocatable, dimension(:,:) :: lhsK real*8, allocatable, dimension(:,:) :: EGmass real*8, allocatable, dimension(:,:) :: EGmasst inquire(file='xyzts.dat',exist=exts) lskeep=lstep if(exts) then open(unit=626,file='xyzts.dat',status='old') read(626,*) ntspts, freq, tolpt, iterat, varcod call sTD ! sets data structures do jj=1,ntspts ! read coordinate data where solution desired read(626,*) ptts(jj,1),ptts(jj,2),ptts(jj,3) enddo close(626) statptts(:,:) = 0 parptts(:,:) = zero varts(:,:) = zero allocate (ivarts(ntspts*ndof)) allocate (ivartsg(ntspts*ndof)) allocate (vartssoln(ntspts*ndof)) allocate (vartssolng(ntspts*ndof)) nbuff=ntout allocate (vartsbuff(ntspts,ndof,nbuff)) if (myrank .eq. master) then do jj=1,ntspts fvarts='varts/varts' fvarts=trim(fvarts)//trim(cname2(jj)) fvarts=trim(fvarts)//trim(cname2(lstep)) fvarts=trim(fvarts)//'.dat' fvarts=trim(fvarts) open(unit=1000+jj, file=fvarts, status='unknown') enddo endif endif c c c.... open history and aerodynamic forces files c if (myrank .eq. master) then open (unit=ihist, file=fhist, status='unknown') open (unit=iforce, file=fforce, status='unknown') endif c c c.... initialize c ifuncs = 0 ! func. evaluation counter istep = 0 ntotGM = 0 ! number of GMRES iterations time = 0 yold = y acold = ac if (mod(impl(1),100)/10 .eq. 1) then c c generate the sparse data fill vectors c allocate (rowp(nshg,nnz)) allocate (colm(nshg+1)) call genadj(colm, rowp, icnt ) ! preprocess the adjacency list nnz_tot=icnt ! this is exactly the number of non-zero ! blocks on this proc allocate (lhsK(nflow*nflow,nnz_tot)) endif if (mod(impl(1),100)/10 .eq. 3) then c c generate the ebe data fill vectors c nedof=nflow*nshape allocate (EGmass(numel,nedof*nedof)) endif cc c.......................................... rerr = zero ybar(:,1:ndof) = y(:,1:ndof) do idof=1,5 ybar(:,ndof+idof) = y(:,idof)*y(:,idof) enddo ybar(:,ndof+6) = y(:,1)*y(:,2) ybar(:,ndof+7) = y(:,1)*y(:,3) ybar(:,ndof+8) = y(:,2)*y(:,3) c......................................... vbc_prof(:,1:3) = BC(:,3:5) c c.... loop through the time sequences c do 3000 itsq = 1, ntseq itseq = itsq c c.... set up the current parameters c nstp = nstep(itseq) nitr = niter(itseq) LCtime = loctim(itseq) c call itrSetup ( y, acold) isclr=0 niter(itseq)=0 ! count number of flow solves in a step ! (# of iterations) do i=1,seqsize if(stepseq(i).eq.0) niter(itseq)=niter(itseq)+1 enddo nitr = niter(itseq) c c.... determine how many scalar equations we are going to need to solve c nsclrsol=nsclr ! total number of scalars solved. At ! some point we probably want to create ! a map, considering stepseq(), to find ! what is actually solved and only ! dimension EGmasst to the appropriate ! size. if(nsclrsol.gt.0)allocate (EGmasst(numel*nshape*nshape & ,nsclrsol)) c c.... loop through the time steps c ttim(1) = REAL(secs(0.0)) / 100. ttim(2) = secs(0.0) c tcorecp1 = REAL(secs(0.0)) / 100. c tcorewc1 = secs(0.0) if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) if(myrank.eq.0) then tcorecp1 = TMRC() endif rmub=datmat(1,2,1) if(rmutarget.gt.0) then rmue=rmutarget xmulfact=(rmue/rmub)**(1.0/nstp) if(myrank.eq.0) then write(*,*) 'viscosity will by multiplied by ', xmulfact write(*,*) 'to bring it from ', rmub,' down to ', rmue endif datmat(1,2,1)=datmat(1,2,1)/xmulfact ! make first step right else rmue=datmat(1,2,1) ! keep constant xmulfact=one endif if(iramp.eq.1) then call initBCprofileScale(vbc_prof,BC,yold,x) ! fix the yold values to the reset BC call itrBC (yold, ac, iBC, BC, iper, ilwork) isclr=1 ! fix scalar call itrBCsclr(yold,ac,iBC,BC,iper,ilwork) endif 867 continue do 2000 istp = 1, nstp if(iramp.eq.1) & call BCprofileScale(vbc_prof,BC,yold) c call rerun_check(stopjob) cc if(stopjob.ne.0) goto 2001 c c Decay of scalars c if(nsclr.gt.0 .and. tdecay.ne.1) then yold(:,6:ndof)=y(:,6:ndof)*tdecay BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*tdecay endif if(nosource.eq.1) BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*0.8 c xi=istp*one/nstp c datmat(1,2,1)=rmub*(1.0-xi)+xi*rmue datmat(1,2,1)=xmulfact*datmat(1,2,1) c.... if we have time varying boundary conditions update the values of BC. c these will be for time step n+1 so use lstep+1 c if(itvn.gt.0) call BCint((lstep+1)*Delt(1), shp, shgl, & shpb, shglb, x, BC, iBC) if(iLES.gt.0) then c c.... get dynamic model coefficient c ilesmod=iLES/10 c c digit bit set filter rule, 10 bit set model c if (ilesmod.eq.0) then ! 0 < iLES < 10 => dyn. model calculated ! at nodes based on discrete filtering call getdmc (yold, shgl, shp, & iper, ilwork, nsons, & ifath, x) endif if (ilesmod .eq. 1) then ! 10 < iLES < 20 => dynamic-mixed ! at nodes based on discrete filtering call bardmc (yold, shgl, shp, & iper, ilwork, & nsons, ifath, x) endif if (ilesmod .eq. 2) then ! 20 < iLES < 30 => dynamic at quad ! pts based on lumped projection filt. call projdmc (yold, shgl, shp, & iper, ilwork, x) endif c endif c c.... set traction BCs for modeled walls c if (itwmod.ne.0) then !wallfn check call asbwmod(yold, acold, x, BC, iBC, & iper, ilwork, ifath, velbar) endif c c.... -----------------------> predictor phase <----------------------- c call itrPredict( yold, acold, y, ac ) call itrBC (y, ac, iBC, BC, iper, ilwork) isclr = zero if (nsclr.gt.zero) then do isclr=1,nsclr call itrBCSclr (y, ac, iBC, BC, iper, ilwork) enddo endif c c.... --------------------> multi-corrector phase <-------------------- c iter=0 ilss=0 ! this is a switch thrown on first solve of LS redistance c cHACK to make it keep solving RANS until tolerance is reached c istop=0 iMoreRANS=0 c c find the the RANS portion of the sequence c do istepc=1,seqsize if(stepseq(istepc).eq.10) iLastRANS=istepc enddo iseqStart= 1 9876 continue c do istepc=iseqStart,seqsize icode=stepseq(istepc) if(mod(icode,10).eq.0) then ! this is a solve isolve=icode/10 if(isolve.eq.0) then ! flow solve (encoded as 0) c etol=epstol(1) iter = iter+1 ifuncs(1) = ifuncs(1) + 1 c c.... reset the aerodynamic forces c Force(1) = zero Force(2) = zero Force(3) = zero HFlux = zero c c.... form the element data and solve the matrix problem c c.... explicit solver c if (impl(itseq) .eq. 0) then if (myrank .eq. master) & call error('itrdrv ','impl ',impl(itseq)) endif if (mod(impl(1),100)/10 .eq. 1) then ! sparse solve c c.... preconditioned sparse matrix GMRES solver c lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) iprec=lhs nedof = nflow*nshape c write(*,*) 'lhs=',lhs call SolGMRs (y, ac, yold, & acold, x, & iBC, BC, & colm, rowp, lhsk, & res, & BDiag, a(mHBrg), a(meBrg), & a(myBrg), a(mRcos), a(mRsin), & iper, ilwork, & shp, shgl, & shpb, shglb, solinc, & rerr) else if (mod(impl(1),100)/10 .eq. 2) then ! mfg solve c c.... preconditioned matrix-free GMRES solver c lhs=0 iprec = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) nedof = 0 call SolMFG (y, ac, yold, & acold, x, & iBC, BC, & res, & BDiag, a(mHBrg), a(meBrg), & a(myBrg), a(mRcos), a(mRsin), & iper, ilwork, & shp, shgl, & shpb, shglb, solinc, & rerr) c else if (mod(impl(1),100)/10 .eq. 3) then ! ebe solve c.... preconditioned ebe matrix GMRES solver c lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) iprec = lhs nedof = nflow*nshape c write(*,*) 'lhs=',lhs call SolGMRe (y, ac, yold, & acold, x, & iBC, BC, & EGmass, res, & BDiag, a(mHBrg), a(meBrg), & a(myBrg), a(mRcos), a(mRsin), & iper, ilwork, & shp, shgl, & shpb, shglb, solinc, & rerr) endif c else ! solve a scalar (encoded at isclr*10) ifuncs(isclr+2) = ifuncs(isclr+2) + 1 etol=epstol(isclr+1) isclr=isolve if((iLSet.eq.2).and.(ilss.eq.0) & .and.(isclr.eq.2)) then ilss=1 ! throw switch (once per step) c c... copy the first scalar at t=n+1 into the second scalar of the c... level sets c y(:,6) = yold(:,6) + (y(:,6)-yold(:,6))/alfi y(:,7) = y(:,6) yold(:,7) = y(:,7) ac(:,7) = zero c call itrBCSclr (y, ac, iBC, BC, iper, ilwork) c c....store the flow alpha, gamma parameter values and assigm them the c....Backward Euler parameters to solve the second levelset scalar c alfit=alfi gamit=gami almit=almi alfi = 1 gami = 1 almi = 1 endif c lhs = 1 - min(1,mod(ifuncs(isclr+2)-1, & LHSupd(isclr+2))) iprec = lhs istop=0 call SolGMRSclr(y, ac, yold, & acold, EGmasst(1,isclr), & x, elDw, & iBC, BC, & rest, & a(mHBrg), a(meBrg), & a(myBrg), a(mRcos), a(mRsin), & iper, ilwork, & shp, shgl, & shpb, shglb, solinc(1,isclr+5)) c endif ! end of scalar type solve c c c.... end of the multi-corrector loop c 1000 continue !check this else ! this is an update (mod did not equal zero) iupdate=icode/10 ! what to update if(iupdate.eq.0) then !update flow call itrCorrect ( y, ac, yold, acold, solinc) call itrBC (y, ac, iBC, BC, iper, ilwork) call tnanq(y, 5, 'y_updbc') c Elaine-SPEBC if((irscale.ge.0).and.(myrank.eq.master)) then call genscale(y, x, iBC) c call itrBC (y, ac, iBC, BC, iper, ilwork) endif else ! update scalar isclr=iupdate !unless if(iupdate.eq.nsclr+1) isclr=0 call itrCorrectSclr ( y, ac, yold, acold, & solinc(1,isclr+5)) if (ilset.eq.2 .and. isclr.eq.2) then fct2=one/almi fct3=one/alfi acold(:,7) = acold(:,7) & + (ac(:,7)-acold(:,7))*fct2 yold(:,7) = yold(:,7) & + (y(:,7)-yold(:,7))*fct3 call itrBCSclr ( yold, acold, iBC, BC, & iper, ilwork) ac(:,7) = acold(:,7)*(one-almi/gami) y(:,7) = yold(:,7) ac(:,7) = zero if (ivconstraint .eq. 1) then & c ... applying the volume constraint c call solvecon (y, x, iBC, BC, & iper, ilwork, shp, shgl) c endif ! end of volume constraint calculations endif call itrBCSclr ( y, ac, iBC, BC, iper, ilwork) endif endif !end of switch between solve or update enddo ! loop over sequence in step if((istop.lt.0).and.(iMoreRANS.lt.5)) then iMoreRANS=iMoreRANS+1 if(myrank.eq.0) write(*,*) 'istop =', istop iseqStart=iLastRANS goto 9876 endif c c Find the solution at the end of the timestep and move it to old c c.... First to reassign the parameters for the original time integrator scheme c if((iLSet.eq.2).and.(ilss.eq.1)) then alfi =alfit gami =gamit almi =almit endif call itrUpdate( yold, acold, y, ac) call itrBC (yold, acold, iBC, BC, iper,ilwork) c Elaine-SPEBC if((irscale.ge.0).and.(myrank.eq.master)) then call genscale(yold, x, iBC) c call itrBC (y, ac, iBC, BC, iper, ilwork) endif do isclr=1,nsclr call itrBCSclr (yold, acold, iBC, BC, iper, ilwork) enddo c istep = istep + 1 c c.... compute boundary fluxes and print out c lstep = lstep + 1 ntoutv=max(ntout,100) if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then c c here is where we save our averaged field. In some cases we want to c write it less frequently if( (mod(lstep, ntoutv) .eq. 0) .and. & ((irscale.ge.0).or.(itwmod.gt.0) .or. & ((nsonmax.eq.1).and.(iLES.gt.0)))) & call rwvelb ('out ', velbar ,ifail) call Bflux (yold, acold, x, & shp, shgl, shpb, & shglb, nodflx, ilwork) endif c c.... end of the NSTEP and NTSEQ loops c c... dump TIME SERIES if (exts) then if (mod(lstep-1,freq).eq.0) then if (numpe > 1) then do jj = 1, ntspts vartssoln((jj-1)*ndof+1:jj*ndof)=varts(jj,:) ivarts=zero enddo do k=1,ndof*ntspts if(vartssoln(k).ne.zero) ivarts(k)=1 enddo call MPI_REDUCE(vartssoln, vartssolng, ndof*ntspts, & MPI_DOUBLE_PRECISION, MPI_SUM, master, & MPI_COMM_WORLD, ierr) call MPI_REDUCE(ivarts, ivartsg, ndof*ntspts, & MPI_INTEGER, MPI_SUM, master, & MPI_COMM_WORLD, ierr) if (myrank.eq.zero) then do jj = 1, ntspts indxvarts = (jj-1)*ndof do k=1,ndof if(ivartsg(indxvarts+k).ne.0) then ! none of the vartssoln(parts) were non zero varts(jj,k)=vartssolng(indxvarts+k)/ & ivartsg(indxvarts+k) endif enddo enddo endif !only on master endif !only if numpe > 1 if( istp.eq. nstp) then !make sure incomplete buffers get purged. icheck=mod(nstp,nbuff) if(icheck.ne.0) nbuff=icheck endif if (myrank.eq.zero) then k=mod(lstep,nbuff) if(k.eq.0) k=nbuff do jj = 1, ntspts vartsbuff(jj,1:5,k)=varts(jj,1:5) enddo if(k.eq. nbuff) then do jj = 1, ntspts ifile = 1000+jj do ibuf=1,nbuff write(ifile,555) lstep-1 -nbuff+ibuf, & (vartsbuff(jj,k,ibuf), k=1,5) ! assuming ndof to be 5 enddo ! buff empty c call flush(ifile) enddo ! jj ntspts endif !only dump when buffer full endif !only on master varts(:,:) = zero ! reset the array for next step 555 format(i6,5(2x,E12.5e2)) endif endif c c.... -------------------> error calculation <----------------- if(ierrcalc.eq.1.or.ioybar.eq.1) then tfact=one/istep do idof=1,ndof ybar(:,idof) =tfact*yold(:,idof) + & (one-tfact)*ybar(:,idof) enddo c....compute average c... ybar(:,ndof+1:ndof+8) is for avg. of square as uu, vv, ww, pp, TT, uv, uw, and vw do idof=1,5 ! avg. of square for 5 flow variables ybar(:,ndof+idof) = tfact*yold(:,idof)**2 + & (one-tfact)*ybar(:,ndof+idof) enddo ybar(:,ndof+6) = tfact*yold(:,1)*yold(:,2) + !uv & (one-tfact)*ybar(:,ndof+6) ybar(:,ndof+7) = tfact*yold(:,1)*yold(:,3) + !uw & (one-tfact)*ybar(:,ndof+7) ybar(:,ndof+8) = tfact*yold(:,2)*yold(:,3) + !vw & (one-tfact)*ybar(:,ndof+8) c... compute err c rerr(:, 7)=rerr(:, 7)+(yold(:,1)-ybar(:,1))**2 c rerr(:, 8)=rerr(:, 8)+(yold(:,2)-ybar(:,2))**2 c rerr(:, 9)=rerr(:, 9)+(yold(:,3)-ybar(:,3))**2 c rerr(:,10)=rerr(:,10)+(yold(:,4)-ybar(:,4))**2 endif c........................................................................... 2000 continue 2001 continue ttim(1) = REAL(secs(0.0)) / 100. - ttim(1) ttim(2) = secs(0.0) - ttim(2) c tcorecp2 = REAL(secs(0.0)) / 100. c tcorewc2 = secs(0.0) if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) if(myrank.eq.0) then tcorecp2 = TMRC() write(6,*) 'T(core) cpu = ',tcorecp2-tcorecp1 endif c call wtime 3000 continue c c.... ----------------------> Post Processing <---------------------- c c.... print out the last step c if ((irs .ge. 1) .and. ((mod(lstep, ntout) .ne. 0) .or. & (nstp .eq. 0))) then if( (mod(lstep, ntoutv) .eq. 0) .and. & ((irscale.ge.0).or.(itwmod.gt.0) .or. & ((nsonmax.eq.1).and.(iLES.gt.0)))) & call rwvelb ('out ', velbar ,ifail) call Bflux (yold, acold, x, & shp, shgl, shpb, & shglb, nodflx, ilwork) endif c if(ierrcalc.eq.1) then c c.....smooth the error indicators c c ! errsmooth is currently available only for incompressible code c do i=1,ierrsmooth c call errsmooth( rerr, x, iper, ilwork, shp, shgl, iBC ) c end do call write_error(myrank, lstep, nshg, 10, rerr ) endif if(ioybar.eq.1) then call write_field(myrank,'a','ybar',4, & ybar,'d',nshg,ndof+8,lstep) endif if(iRANS.lt.0 .and. idistcalc.eq.1) then call write_field(myrank,'a','DESd',4, & elDw,'d',numel,1,lstep) call write_field(myrank,'a','dwal',4, & d2wall,'d',nshg,1,lstep) endif c c.... close history and aerodynamic forces files c if (myrank .eq. master) then close (ihist) close (iforce) if(exts) then deallocate(ivarts) deallocate(ivartsg) deallocate(vartssoln) deallocate(vartssolng) do jj=1,ntspts close(1000+jj) enddo endif endif close (iecho) if(iabc==1) deallocate(acs) c c.... end c return end