C----------------------------------------------------------------------- C nek5000 user-file template C C user specified routines: C - userbc : boundary conditions C - useric : initial conditions C - uservp : variable properties C - userf : local acceleration term for fluid C - userq : local source term for scalars C - userchk: 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.0 utrans = 0.0 return end c----------------------------------------------------------------------- subroutine userf (ix,iy,iz,eg) c c Note: this is an acceleration term, NOT a force! c Thus, ffx will subsequently be multiplied by rho(x,t). c include 'SIZE' include 'TOTAL' include 'NEKUSE' integer e,f,eg c e = gllel(eg) 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' call walldist(igeom) !Fills PS1 with wall distance and PS2 with Y+ c-----File IO ifxyo = .true. !Turn on geometry & scalar output ifpsco(1) = .true. ifpsco(2) = .true. ifpsco(3) = .true. call prepost(.true.,'tes') !File prefix 'avg' return end c----------------------------------------------------------------------- subroutine userbc (ix,iy,iz,iside,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' ux = -(y-1.5)*(y-1.5)+1.5*1.5 uy = 0.0 uz = 0.0 temp = 0.0 return end c----------------------------------------------------------------------- subroutine useric (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' ux = -(y-1.5)*(y-1.5)+1.5*1.5 uy = 0.0 uz = 0.0 temp = 0.0 return end c----------------------------------------------------------------------- subroutine usrdat ! This routine to modify element vertices include 'SIZE' include 'TOTAL' return end c----------------------------------------------------------------------- subroutine usrdat2 ! This routine to modify mesh coordinates include 'SIZE' include 'TOTAL' return end c----------------------------------------------------------------------- subroutine usrdat3 include 'SIZE' include 'TOTAL' return end c----------------------------------------------------------------------- subroutine walldist (igeom) !This is more or less from TWALLUZ in turb.f, which in turn is called by drive2.f C include 'SIZE' include 'TOTAL' parameter (lr=lx1*ly1*lz1) common /SCRMG/ XWLL(LX1,LZ1) $ , YWLL(LX1,LZ1) $ , ZWLL(LX1,LZ1) common /SCRCH/ UTW(LX1),ZNW(LX1),USTR(LX1),ZSTR(LX1) real wdist (lx1,ly1,lz1,lelv) $ ,strsinel(lx1,ly1,lz1,lelv) integer iel,ifc character cb*3 c if (igeom.eq.1) return c nface=2*ndim ntot1=nx1*ny1*nz1*nelv call rzero(wdist,ntot1) do iel=1,nelv do ifc=1,nface cb = cbc(ifc,iel,1) IF (cb.eq.'W ') then !Changed BC by Markus call facexv (XWLL,YWLL,ZWLL,XM1(1,1,1,IEL),YM1(1,1,1,IEL), $ ZM1(1,1,1,IEL),IFC,0) call helpwuz (XWLL,YWLL,ZWLL,IEL,IFC) call getstrselem(strsinel(1,1,1,iel),ifc,iel) !strsinel is only nonzero on wall faces !returns sqrt(tau_wall) call facev(wdist,iel,ifc,znw(2),lx1,ly1,lz1) !Fills variable wdist with wall distance endif enddo enddo call copy(t(1,1,1,1,2),wdist,ntot1) !Copy wall distance to PS1 !wdist is constant on each element face call cmult(strsinel,1/(param(1)**0.5*param(2)),ntot1) ! sqrt(tau_wall)/[sqrt(rho)*visc)] call col2(strsinel,wdist,ntot1) ! sqrt(tau_wall)/[sqrt(rho)*visc)]*y call copy(t(1,1,1,1,3),strsinel,ntot1) !Copy y+ to PS2 return end c----------------------------------------------------------------------- SUBROUTINE HELPWUZ (XWLL,YWLL,ZWLL,IEL,IFC) !This is from COMWUZ in turb.f C INCLUDE 'SIZE' INCLUDE 'INPUT' INCLUDE 'SOLN' INCLUDE 'TSTEP' INCLUDE 'GEOM' INCLUDE 'TOPOL' COMMON /SCRCH/ UTW(LX1),ZNW(LX1),USTR(LX1),ZSTR(LX1) C DIMENSION XWLL(LX1,LZ1),YWLL(LX1,LZ1),ZWLL(LX1,LZ1) $ , ISHDAT(6) C DATA ISHDAT / 1,-1,-1, 1, 1,-1/ C NXM1 = NX1 - 1 NXM2 = NX1 - 2 ISCH = ISHDAT(IFC) C CALL FACIND2 (JS1,JF1,JSKIP1,JS2,JF2,JSKIP2,IFC) C DO 500 JWY=1,NZ1 DO 500 JWX=1,NX1 IWX = JS1 + (JWX-1)*JSKIP1 IWY = JS2 + (JWY-1)*JSKIP2 CALL NORIND (JS3,JF3,JSKIP3,IWX,IWY,IFC,ISCH) CALL CWREF (XWLL,YWLL,ZWLL,UTW,ZNW,0,0,0,IEL,IFC, $ JWX,JWY,JS3,JSKIP3) C 500 CONTINUE C RETURN END c----------------------------------------------------------------------- subroutine getstrselem(strsel,f,e) !This computes the square root of wall shear stress !on each grid point of a given face c include 'SIZE' include 'TOTAL' c real strsel(lx1,ly1,lz1) $ ,visc(lx1,ly1,lz1,lelv) $ ,gije(lx1,ly1,lz1,3,3) $ ,sn(3) $ ,StrX,StrY,StrZ integer i,j,k,e,f $ ,kx1,kx2,ky1,ky2,kz1,kz2 nxyz=lx1*ly1*lz1 !Zero out strsel call rzero(strsel,lx1*ly1*lz1) call comp_gije(gije,vx(1,1,1,e),vy(1,1,1,e),vz(1,1,1,e),e) !Compute gradient tensor for entire element call cmult(gije,param(2),lx1*ly1*lz1*9) !Multiply gradient with viscosity if (if3d) then call facind(kx1,kx2,ky1,ky2,kz1,kz2,nx1,ny1,nz1,f) do i=kx1,kx2 do j=ky1,ky2 do k=kz1,kz2 call getSnormal(sn,i,j,k,f,e) !Compute surface normal vector StrX=gije(i,j,k,1,1)*sn(1)+gije(i,j,k,1,2)*sn(2) !Projection and summation & +gije(i,j,k,1,3)*sn(3) StrY=gije(i,j,k,2,1)*sn(1)+gije(i,j,k,2,2)*sn(2) & +gije(i,j,k,2,3)*sn(3) StrZ=gije(i,j,k,3,1)*sn(1)+gije(i,j,k,3,2)*sn(2) & +gije(i,j,k,3,3)*sn(3) strsel(i,j,k)=(StrX**2+StrY**2+StrZ**2)**0.25 !Note: this computes sqrt(tau_wall) enddo enddo enddo else print *,'Code not made for 2D/Axi' call exitt endif 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