c----------------------------------------------------------------------- c This is same as gs2a, save the geometry is modified to accomodate c a smooth transition at the channel inlet. c c The transition is realized by the following actions: c c 1. Rescale in y to expand the geometry to +/- (1+h_in) c c 2. Shift geometry by -1 unit in x, so nozzle region will be on [-1,0] c c 3. Shape nozzle region to transition from +/-(1+h_in) to +/-1 on [-1,0] c c 4. Rescale y for all x > 0 to be on +/-1. c c 5. Develop grooves as before, but now on xm1,ym1 in usrdat2, c instead of on xc,yc in usrdat. c c Note: For the grooved channel formation that is generated in usrdat c or usrdat2, we assume the following about the original geometry: c c c 1. Generated by genbox c c 2. Test section is in x > 0 c c 3. There are an even number of elements on each interval c c [0,2] [2,4] [4,6] ... [xn-2,xn] c c and that such a region, on [0,xn] will in the end yield n/2 grooves c c 4. If you want more grooves, increase xn in the gn2a.box file that is c the input to "genbox" c c 5. If you want more resolution in the x-direction, increase "nelx" in c gn2a.box for Box 1, which is the test section (plus smooth transition) 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' real lambda 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 c ubar = 1. c call sponge(x,lambda) c ffx = lambda*(ubar-ux) ! sponge 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 return end c----------------------------------------------------------------------- subroutine userchk include 'SIZE' include 'TOTAL' common /scruz/ wt(lx1*ly1*lz1,lelv) integer e nxyz = nx1*ny1*nz1 ! # pts in an element call rzero(wt,n) call fill_div(usrdiv) ! for outflow conditions xmin = 1.e33 xmax = -1.e33 do e=1,nelv x = xm1(2,2,1,e) if (0.lt.x.and.x.lt.2) then xmnl = vlmin(xm1(1,1,1,e),nxyz) xmxl = vlmax(xm1(1,1,1,e),nxyz) xmin = min(xmin,xmnl) xmax = max(xmax,xmxl) call rone(wt(1,e),nxyz) ! add this el to list endif enddo xmin = glmin(xmin,1) xmax = glmax(xmax,1) ! ubar = glsc3(vx,bm1,wt,n)/glsc2(bm1,wt,n) ! Conditional average vol_nominal = 2*(xmax-xmin) ! Nominal channel volume, assuming ubar = glsc3(vx,bm1,wt,n)/vol_nominal ! channel 1/2-height=1 c param(68) = 100 c call avg_all n = nx1*ny1*nz1*nelt tmin = glmin(t ,n) tmax = glmax(t ,n) umin = glmin(vx,n) umax = glmax(vx,n) if (mod(istep,10).eq.0.and.nid.eq.0) $ write(6,1) istep,time,umin,umax,tmin,tmax 1 format(i9,1p5e13.5,' tmax') return end c----------------------------------------------------------------------- subroutine userbc (ix,iy,iz,iside,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' common /mycoords/ shiftx xx = x-shiftx ! Effective origin of inflow r2 = xx*xx+y*y scale = 4./pi ux = -scale*xx/r2 ! Radial inflow from far field uy = -scale*y/r2 uz = 0.0 temp=0.0 if (x.gt.0) temp=10.0 return end c----------------------------------------------------------------------- subroutine useric (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' common /mycoords/ shiftx ux=0.0 uy=0.0 uz=0.0 temp=0 if (x.lt.shiftx) then xx = x-shiftx ! Effective origin of inflow r2 = xx*xx+y*y scale = 2./pi if (r2.gt.0) then ux = -scale*xx/r2 ! Radial inflow from far field uy = -scale*y/r2 endif else ux = 1 endif return end c----------------------------------------------------------------------- subroutine usrdat2 include 'SIZE' include 'TOTAL' parameter (mgrooves=200) real gdx(0:mgrooves),xg(0:mgrooves),ght(0:mgrooves) n = nx1*ny1*nz1*nelt c This is the initial geometry translation to set up for smooth inlet h_in = 0.75 scaley = 1+h_in call cmult(ym1,scaley,n) shiftx = -1 call cadd (xm1,shiftx,n) xmin = glmin(xm1,n) xmax = glmax(xm1,n) if (nid.eq.0) write(6,*) xmin,xmax,' xmax0' c This is the nozzle transition. c We choose it to be a sine, but anything reasonable is ok. do i=1,n x = xm1(i,1,1,1) y = ym1(i,1,1,1) if (x.gt.0) then ! in the test section ynew = y / scaley elseif (x.gt.shiftx) then argx = pi*(x-shiftx)/(2*abs(shiftx)) disp = -h_in*sin(argx) * (y/scaley) ynew = y+disp else ynew = y endif ym1(i,1,1,1) = ynew enddo c Set up groove parameters and modify test section: x>0 ngrooves = 80 ! Assumes 40 elements of length 1 in x-dir, x>0 ght (0) = 1.20 ! Initial groove height gdx (0) = 0.60 ! Initial groove size xg (0) = 0 do i=1,ngrooves ig1 = min(40,i-1) xg (i) = xg (i-1) + gdx(i-1) ght(i) = ght(i-1) gdx(i) = 1.03*gdx(ig1) enddo do i=1,n x = xm1(i,1,1,1) y = ym1(i,1,1,1) if (x.gt.0) then ! in the test section ig = x/2. ! Map original coord. to groove number ss = (x-2*ig)/2. if (ss.le.0.5) then disp = ght(ig)*y*ss ! Ascending part of groove else disp = ght(ig)*y*(1-ss) ! Descending part of groove endif xnew = xg(ig) + ss*gdx(ig) ynew = y + disp xm1(i,1,1,1) = xnew ym1(i,1,1,1) = ynew c write(6,1) i,ig,x,y,xnew,ynew c 1 format(2i8,1p4e14.6,' xnew') endif enddo param(59) = 1. ! Set all elements to deformed return end c----------------------------------------------------------------------- subroutine usrdat include 'SIZE' include 'TOTAL' param(66) = 4. ! These give the std nek binary i/o and are param(67) = 4. ! good default values return end c----------------------------------------------------------------------- subroutine usrdat3 include 'SIZE' include 'TOTAL' c return end c----------------------------------------------------------------------- subroutine sponge(x,lambda) include 'SIZE' real x,xstart,xend,delta_rise,delta_fall real lambda,lambda_max,s1,s2,arg1,arg2 integer icalld save icalld data icalld /0/ real xmax save xmax if (icalld.eq.0) then icalld=1 n=nx1*ny1*nz1*nelv xmax = glmax(x,n) endif xstart = 0.975*xmax xend = xmax delta_rise = 1.5 delta_fall = 0.5 lambda_max = 1. arg1 = (x-xstart)/delta_rise arg2 = ((x-xend)/delta_fall) + 1 call smooth(s1,arg1) call smooth(s2,arg2) lambda = lambda_max*(s1-s2) return end c----------------------------------------------------------------------- subroutine smooth(s,x) real s,x if (x .lt. 0.) then s = 0. elseif (x .gt. 0. .and. x .lt. 1.) then s = 1./(1.+exp(1./(x-1.)+1./x)) elseif (x .ge. 1) then s = 1. endif 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' common /cvflow_d/ dist(lx2,ly2,lz2,lelt),dmax integer e,f nxyz2 = nx2*ny2*nz2 ntot2 = nx2*ny2*nz2*nelv ntot1 = nx1*ny1*nz1*nelv call rzero(dist,ntot2) davg = 0. wavg = 0. xmax = glmax(xm1,ntot1) nface = 2*ndim do e=1,nelv do f=1,nface if (cbc(f,e,1).eq.'O ') then ! Outflow do k=1,nxyz2 dist(k,1,1,e) = xmax-xm2(k,1,1,e) davg = davg+dist(k,1,1,e) wavg = wavg+1.0 enddo endif enddo 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-----------------------------------------------------------------------