c----------------------------------------------------------------------- C C USER SPECIFIED ROUTINES: C C - boundary conditions C - initial conditions C - variable properties C - local acceleration for fluid (a) C - forcing function for passive scalar (q) C - general purpose routine for checking errors etc. C c----------------------------------------------------------------------- subroutine uservp (ix,iy,iz,eg) include 'SIZE' include 'TOTAL' include 'NEKUSE' integer e,f,eg c e = gllel(eg) udiff =0. utrans=0. return end c----------------------------------------------------------------------- subroutine userf (ix,iy,iz,eg) include 'SIZE' include 'TOTAL' include 'NEKUSE' integer e,f,eg c e = gllel(eg) c Note: this is an acceleration term, NOT a force! c Thus, ffx will subsequently be multiplied by rho(x,t). ffx = 0.0 ffy = 0.0 ffz = 0.0 return end c----------------------------------------------------------------------- subroutine userq (ix,iy,iz,eg) include 'SIZE' include 'TOTAL' include 'NEKUSE' integer e,f,eg c e = gllel(eg) qvol = 0.0 source = 0.0 return end c----------------------------------------------------------------------- subroutine userchk include 'SIZE' include 'TOTAL' real linf,l2,h1,semi,Energy c call hpts() call fill_div(usrdiv) if (istep.gt.0) then call normvc(h1,semi,l2,linf,vx,vy,vz) Energy = .5*l2*l2 if (nid.eq.0) write(06,*) istep,time,Energy,' energy' endif if (istep.gt.0.and.istep.lt.nbdinp) call my_full_restart_load call my_full_restart_save ! save add'l files for full-restart return end c----------------------------------------------------------------------- subroutine userbc (ix,iy,iz,iside,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' rr = x*x + y*y ux = 0.0 uy = 0.0 uz = 2*(1 - 4*rr) ! Hagen-Poiseiulle flow, pipe diamter = 1. temp= 0.0 return end c----------------------------------------------------------------------- subroutine useric (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' rr = x*x +y*y ux=0.0 uy=0.0 uz=2*(1 - 4*rr) temp=0 return end c----------------------------------------------------------------------- subroutine usrdat include 'SIZE' include 'TOTAL' c return end c----------------------------------------------------------------------- subroutine usrdat3 include 'SIZE' include 'TOTAL' c return end c----------------------------------------------------------------------- subroutine usrdat2 ! Modify geometry include 'SIZE' include 'TOTAL' logical ifstenosis save ifstenosis c param(66) = 4. c param(67) = 4. x0 = -.5 x1 = .5 call rescale_x (xm1,x0,x1) ! Rescale x & y so pipe diameter is 1. call rescale_x (ym1,x0,x1) ifstenosis = .false. ifstenosis = .true. if (.not.ifstenosis) return one = 1.0 pi = 4.*atan(one) n = nx1*ny1*nz1*nelv z_outflow = 45. z_inflow = -5.0 call rescale_x(zm1,z_inflow,z_outflow) stenosis_length = 1. do i=1,n x = xm1(i,1,1,1) y = ym1(i,1,1,1) z = zm1(i,1,1,1)/stenosis_length arg = pi amp = .25 ! 75% stenosis for 3D case if (abs(z).le.1) then xm1(i,1,1,1) = x*(1.-amp*(1.+cos(arg*z))) ym1(i,1,1,1) = y*(1.-amp*(1.+cos(arg*z))) c write(6,1) i,x,y,z,xm1(i,1,1,1),ym1(i,1,1,1),zm1(i,1,1,1) endif enddo param(59) = 1 ! Force ifdfrm=.true. ( 8/26/03 ) ifxyo = .true. c call outpost(xm1,ym1,zm1,pr,t,' ') c call exitt return end c----------------------------------------------------------------------- subroutine fill_div(div) c c Fill the domain with a nontrivial divergence, where desired. c c include 'SIZE' include 'TOTAL' c integer icalld save icalld data icalld /0/ c common /cvflow_d/ dist(lx2,ly2,lz2,lelt),dmax common /cvflow_l/ ifdivf logical ifdivf c real div(lx2,ly2,lz2,lelv) c ntot2 = nx2*ny2*nz2*nelv call rzero(div,ntot2) if (ifdivf) return c if (icalld.eq.0) then icalld = 1 call set_outflow_dist endif dmax2 = dmax*dmax c call get_div_const(cdiv) c nd = 0 do i=1,ntot2 dd = dist(i,1,1,1)*dist(i,1,1,1) dd = min(dmax2,dd) if (dist(i,1,1,1).ne.0) div(i,1,1,1) = cdiv*(1-dd/dmax2) if (dist(i,1,1,1).ne.0) nd = nd+1 enddo dsmin = glmin(dist,ntot2) dsmax = glmax(dist,ntot2) dvmin = glmin(div,ntot2) dvmax = glmax(div,ntot2) xd = nd xd = glsum(xd,1) nd = xd if(nid.eq.0)write(6,1) istep,nd,dvmin,dvmax,dsmin,dsmax,cdiv,dmax2 1 format(2i9,1p6e11.3,'divmnmx') 2 format(2i9,1p2e11.3,'divdstm') c return end c----------------------------------------------------------------------- subroutine set_outflow_dist c c Compute projected normal distance from outflow c include 'SIZE' include 'TOTAL' c common /cvflow_d/ dist(lx2,ly2,lz2,lelt),dmax c integer e,f c c nxyz2 = nx2*ny2*nz2 ntot2 = nx2*ny2*nz2*nelv ntot1 = nx1*ny1*nz1*nelv c call rzero(dist,ntot2) davg = 0. wavg = 0. c zmax = glmax(zm1,ntot1) c do e=1,nelv if (cbc(6,e,1).eq.'O ') then ! Outflow do k=1,nxyz2 dist(k,1,1,e) = zmax-zm2(k,1,1,e) davg = davg+dist(k,1,1,e) wavg = wavg+1.0 enddo endif enddo dmax = glmax(dist,ntot2) davg = glsum(davg,1) wavg = glsum(wavg,1) if (wavg.gt.0) davg = davg/wavg if (nid.eq.0) write(6,1) dmax,davg,wavg 1 format('div: davg:',1p3e12.4) dmax = 0.5*(dmax+davg) c return end c----------------------------------------------------------------------- subroutine get_div_const(cdiv) c c Get constant multiplier for divergence c include 'SIZE' include 'TOTAL' c common /cvflow_d/ dist(lx2,ly2,lz2,lelt),dmax real flx(0:2) c integer e,f vnmo_desired = 2. cdivt = 1.5*vnmo_desired/dmax cdiv = max(0.,cdivt) ! No contraction! if (nid.eq.0) write(6,1) istep,time,vnmo,dmax,cdiv,cdivt 1 format(i9,1p5e13.5,' cdiv') return end c----------------------------------------------------------------------- subroutine my_full_restart_save c Saves files for next full restart include 'SIZE' include 'TOTAL' include 'RESTART' ! max_rst c This is the full-restart save part: max_rst = 2*(nbdinp-1) ! max # of rst files saved nps1 = 0 if (ifheat) nps1 = 1 + npscal mostep = mod1(istep,iostep) if (istep.gt.iostep.and.mostep.lt.nbdinp) $ call outpost2(vx,vy,vz,pr,t,nps1,'rst') return end c----------------------------------------------------------------------- subroutine my_full_restart_load include 'SIZE' include 'TOTAL' character*80 s80 call blank(s80,80) if (istep.eq.1) s80 ='rststenosis0.f00001' ! This would be the case for Torder=3 if (istep.eq.2) s80 ='rststenosis0.f00002' ! using the first restart pair c if (istep.eq.1) s80 ='rststenosis.fld03' ! This would be the case for Torder=3 c if (istep.eq.2) s80 ='rststenosis.fld04' ! using the second restart pair c if (istep.eq.1) s80 ='rststenosis0.f00001' ! This would be the case for Torder=2, file 1 c if (istep.eq.1) s80 ='rststenosis0.f00002' ! This would be the case for Torder=2, file 2 call bcast(s80,80) call chcopy(initc,s80,80) c time_curr = time nfiles = 1 call restart(nfiles) ! Note -- time is reset. if (nid.ne.0) time=0 time = glmax(time,1) ! Synchronize time across all processors c time = time_curr ! Preserve current simulation time return end c----------------------------------------------------------------------- c automatically added by makenek subroutine usrsetvert(glo_num,nel,nx,ny,nz) ! to modify glo_num integer*8 glo_num(1) return end