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' include 'ZPER' integer e,f c---------------------------------------------------------------| c ** Routine to calculate the gradients - 05/29/10 ** c---------------------------------------------------------------| parameter(lt=lx1*ly1*lz1) common /scrns/ pm1(lx1,ly1,lz1,lelv) $,vxx(lx1,ly1,lz1,lelv),vxy(lx1,ly1,lz1,lelv),vxz(lx1,ly1,lz1,lelv) $,vyx(lx1,ly1,lz1,lelv),vyy(lx1,ly1,lz1,lelv),vyz(lx1,ly1,lz1,lelv) common /scruz/ $ vzx(lx1,ly1,lz1,lelv),vzy(lx1,ly1,lz1,lelv),vzz(lx1,ly1,lz1,lelv) $,temp(4,ly1**2),work(4,ly1**2) if (istep.eq.0) call set_obj ! obj for surface integrals. if (nid.eq.0.and. istep.eq.0) then ! ONLY NODE 0 opens file open(unit=103,file='coord.txt') open(unit=104,file='elenos.txt') endif call gradm1(vxx,vxy,vxz,vx) ! MORE EFFICIENT than dudxyz DO 100 II=1,NHIS IF (HCODE(10,II).NE.'I') GOTO 100 IOBJ = LOCHIS(1,II) MEMTOT = NMEMBER(IOBJ) IF (HCODE(1,II).NE.' ' .OR. HCODE(2,II).NE.' ' .OR. $ HCODE(3,II).NE.' ' ) THEN IFIELD = 1 DO 50 MEM=1,MEMTOT ISK = 0 IEG = OBJECT(IOBJ,MEM,1) IFC = OBJECT(IOBJ,MEM,2) if (nid.eq.0) write(104,*) ieg,gllel(ieg),gllnid(ieg) call rzero(temp,4*ly1*ly1) IF (gllnid(ieg).eq.nid) then ! this processor has a contribution ie = gllel(ieg) CALL DSSET(NX1,NY1,NZ1) IFACE = EFACE1(IFC) JS1 = SKPDAT(1,IFACE) JF1 = SKPDAT(2,IFACE) JSKIP1 = SKPDAT(3,IFACE) JS2 = SKPDAT(4,IFACE) JF2 = SKPDAT(5,IFACE) JSKIP2 = SKPDAT(6,IFACE) I = 0 DO 110 J2=JS2,JF2,JSKIP2 DO 110 J1=JS1,JF1,JSKIP1 I = I+1 temp(1,i) = xm1(j1,j2,1,ie) temp(2,i) = ym1(j1,j2,1,ie) temp(3,i) = zm1(j1,j2,1,ie) temp(4,i) = vxy(j1,j2,1,ie) 110 CONTINUE ENDIF call gop(temp,work,'+ ',4*ly1*ly1) ! gather results onto rank 0 IF (nid.eq.0) write(103,18) ((temp(k,i),k=1,4),i=1,ly1*ly1) 18 format(1p4e15.7) 50 CONTINUE ENDIF 100 CONTINUE c--------------------------------------------------------------------------- c----------------------------------------------------------| c *** Adding Divergence at the outflow *** c----------------------------------------------------------| nxyz = nx1*ny1*nz1 n = nxyz*nelv if (istep.gt.iostep) ifxyo = .false. ! Turn off xyz output umx = glmax(vx,n) vmx = glmax(vy,n) uix = -9.e20 ! Find max part of domain where div=0 vix = -9.e20 nface = 2*ndim do e=1,nelv iflag=1 do f=1,nface if (cbc(f,e,1).eq.'O ') iflag=0 ! outflow element enddo if (iflag.eq.1) then uit=vlmax(vx(1,1,1,e),nxyz) vit=vlmax(vy(1,1,1,e),nxyz) uix= max(uix,uit) vix= max(vix,vit) endif enddo uix = glmax(uix,1) vix = glmax(vix,1) if (nid.eq.0) write(6,1) istep,time,umx,vmx,uix,vix 1 format(i8,1p5e12.4,' umx') call fill_div(usrdiv) ! nontrivial divergence at outflow return end c----------------------------------------------------------------------- subroutine userbc (ix,iy,iz,iside,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' real xx,yy,lenz,leny,scalerms one = 1.0 pi = 4.*atan(one) lenz = 0.75 ! JS : From Visit. ( Length of the periodic box ) leny = 1.0 scalerms = 0.046 ! such that u_rms/U_infty=0.25, from Matlab script xx = z*2*pi/lenz yy = y*2*pi/leny u = cos(3.0*xx)*cos(4.0*yy)*2*pi/lenz+sin(5.0*yy)*2*pi/lenz v = 0.75*sin(3.0*xx)*sin(4.0*yy)*2*pi/leny+cos(5.0*xx)*2*pi/leny c--------------------------------------------------- c uy = -u*cos(55) ; ux = u*sin(55) | c--------------------------------------------------- ux = 0.81915 uy = v*scalerms*2 - 0.5735 uz = u*scalerms*2 temp = 0 return end c----------------------------------------------------------------------- subroutine useric (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' ux = 0.81915 uy = -0.5735 uz = 0.0 temp = 0 return end c----------------------------------------------------------------------- subroutine usrdat include 'SIZE' include 'TOTAL' return end c----------------------------------------------------------------------- subroutine usrdat2 include 'SIZE' include 'TOTAL' ntot = nx1*ny1*nz1*nelv fac = 0.6 call cmult(zm1, fac, ntot) return end c----------------------------------------------------------------------- subroutine usrdat3 include 'SIZE' include 'TOTAL' c return end c----------------------------------------------------------------------- subroutine fill_div(div) c Fill the domain with a nontrivial divergence, where desired. 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 c vnmo_desired = 3. vnmo_desired = 1.5 cdivt = 1.5*vnmo_desired/dmax cdiv = max(0.,cdivt) ! No contraction! if (nid.eq.0) write(6,1) istep,time,vnmo_desired,dmax,cdiv,cdivt 1 format(i9,1p5e13.5,' cdiv') return end c----------------------------------------------------------------------- subroutine set_obj ! define objects for surface integrals include 'SIZE' include 'TOTAL' integer e,f c Define new objects nobj = 1 iobj = 0 do ii=nhis+1,nhis+nobj iobj = iobj+1 hcode(10,ii) = 'I' hcode( 1,ii) = 'F' ! 'F' hcode( 2,ii) = 'F' ! 'F' hcode( 3,ii) = 'F' ! 'F' lochis(1,ii) = iobj enddo nhis = nhis + nobj if (maxobj.lt.nobj) write(6,*) 'increase maxobj in SIZEu. rm *.o' if (maxobj.lt.nobj) call exitt c nxyz = nx1*ny1*nz1 p = -0.25 do e=1,nelv do f=1,2*ndim if (cbc(f,e,1).eq.'W '.and. xm1(1,1,1,e).gt.p ) then ieg = lglel(e) iobj = 0 if (f.ge.2 .and. f.le.4 ) iobj=1 ! lower wall c Faces vary bet 2 and 4 as found from usrchk subroutine. if (iobj.gt.0) then c write(10,*)f,e nmember(iobj) = nmember(iobj) + 1 mem = nmember(iobj) ieg = lglel(e) object(iobj,mem,1) = ieg object(iobj,mem,2) = f c write(6,1) iobj,mem,f,ieg,e,nid,' OBJ' c write(10,*) iobj,mem,f,ieg,e,nid,'OBJ' 1 format(6i9,a4) endif c endif enddo enddo c write(6,*) 'number',(nmember(k),k=1,4) return end c--------------------------------------------------- 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