c----------------------------------------------------------------------- subroutine uservp(i,j,k,eg) ! set variable properties include 'SIZE' include 'TOTAL' include 'NEKUSE' integer i,j,k,e,eg c e = gllel(eg) udiff = 0.0 utrans = 0.0 return end c----------------------------------------------------------------------- subroutine userf(i,j,k,eg) ! set acceleration term 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 i,j,k,e,eg c e = gllel(eg) ffx = 0.0 ffy = 0.0 ffz = 0.0 return end c----------------------------------------------------------------------- subroutine userq(i,j,k,eg) ! set source term include 'SIZE' include 'TOTAL' include 'NEKUSE' integer i,j,k,e,eg c e = gllel(eg) qvol = 0.0 return end c----------------------------------------------------------------------- subroutine userbc(i,j,k,f,eg) ! set up boundary conditions c NOTE ::: This routine MAY NOT be called by every process include 'SIZE' include 'TOTAL' include 'NEKUSE' integer i,j,k,f,e,eg e=gllel(eg) ux = 0.0 uy = 0.0 uz = 1.0 if (cbu.eq.'o ') then U0 = 1.0 ! characteristic velocity delta = 0.1 ! small positive constant pa = dongOutflow(i,j,k,e,f,U0,delta) endif flux = 0.0 temp = 0.0 c if (cbc(f,e,1).eq.'W '.AND.cbc(f,e,2).eq.'t ') temp=1.0 if (cbc(f,e,1).eq.'W '.AND.cbc(f,e,2).eq.'f ') flux=1.0 return end c----------------------------------------------------------------------- subroutine useric(i,j,k,eg) ! set up initial conditions include 'SIZE' include 'TOTAL' include 'NEKUSE' integer i,j,k,e,eg ux = 0.0 uy = 0.0 uz = 1.0 temp = 0.0 return end c----------------------------------------------------------------------- subroutine userchk include 'SIZE' include 'TOTAL' n=lx1*ly1*lz1*nelt if (istep.eq.0) then xxmax = glmax(xm1,n) yymax = glmax(ym1,n) zzmax = glmax(zm1,n) xxmin = glmin(xm1,n) yymin = glmin(ym1,n) zzmin = glmin(zm1,n) if (nio.eq.0) then write(6,18) xxmin,yymin,zzmin write(6,19) xxmax,yymax,zzmax 18 format(' xyz min ',5g13.5) 19 format(' xyz max ',5g13.5) endif ifxyo=.true. c call outpost(vx,vy,vz,t,pr,' ') endif c if (mod(istep,500).eq.0) then c umin=glmin(vx,n) c umax=glmax(vx,n) c vmin=glmin(vy,n) c vmax=glmax(vy,n) c wmin=glmin(vz,n) c wmax=glmax(vz,n) c tmin=glmin(t ,n) c tmax=glmax(t ,n) c if (nio.eq.0) write(6,1) c \$ istep,time,umin,umax,vmin,vmax,wmin,wmax,tmin,tmax c 1 format(i9,1p9e11.3,' tmax') c c call compute_cfl(cfl,vx,vy,vz,dt) c call outpost(cflf,vy,vz,pr,t,'cfl') c endif c if (istep.gt.0) then c ifield = 1 ! essential for lambda2 right now c call lambda2(t(1,1,1,1,2)) c tmin=glmin(t(1,1,1,1,2) ,n) c tmax=glmax(t(1,1,1,1,2) ,n) c if (nio.eq.0) write(6,2)istep,time,tmin,tmax c 2 format(i9,1p3e11.3,' tmax2') c c ifxyo = .true. c ifvo = .true. c ifpro = .true. c ifto = .true. c ifpsco(1) = .true. cc call outpost(vx,vy,vz,pr,t,' ') c endif c call comp_Reynolds c if (mod(istep,500).ne.0) then c call compute_cfl(cfl,vx,vy,vz,1.0) c call outpost(cflf,vy,vz,pr,t,'cfl') c endif return end c----------------------------------------------------------------------- subroutine usrdat ! This routine to modify element vertices include 'SIZE' include 'TOTAL' return end c----------------------------------------------------------------------- subroutine usrdat3 include 'SIZE' include 'TOTAL' return end c----------------------------------------------------------------------- subroutine usrdat2 include 'SIZE' include 'TOTAL' integer e,f ifheat = .false. nface = 2*ldim do e=1,nelv do f=1,nface if(cbc(f,e,1).eq.'E '.or.cbc(f,e,1).eq.' ')then cbc(f,e,1)='E ' cbc(f,e,2)='E ' elseif(cbc(f,e,1).eq.'W ')then ! pebbles boundaryID(f,e) = 4 cbc(f,e,1)='W ' cbc(f,e,2)='f ' elseif(cbc(f,e,1).eq.'W01')then ! side of cylinder inner boundaryID(f,e) = 3 cbc(f,e,1)='W ' cbc(f,e,2)='I ' elseif(cbc(f,e,1).eq.'W02')then ! side of cylinder outer boundaryID(f,e) = 3 cbc(f,e,1)='W ' cbc(f,e,2)='I ' elseif(cbc(f,e,1).eq.'W03')then ! z_bottom boundaryID(f,e) = 1 cbc(f,e,1)='v ' ! Dong's outflow cbc(f,e,2)='t ' elseif(cbc(f,e,1).eq.'W04')then ! z_top boundaryID(f,e) = 2 cbc(f,e,1)='o ' ! Dong's outflow cbc(f,e,2)='O ' else write(*,*)'bc incorrect',e,f,cbc(f,e,1) endif enddo enddo param(18) = 1 ! use algebraic residual norm return end c----------------------------------------------------------------------- function dongOutflow(ix,iy,iz,iel,iside,u0,delta) include 'SIZE' include 'SOLN' include 'GEOM' real sn(3) ux = vx(ix,iy,iz,iel) uy = vy(ix,iy,iz,iel) uz = vz(ix,iy,iz,iel) call getSnormal(sn,ix,iy,iz,iside,iel) vn = ux*sn(1) + uy*sn(2) + uz*sn(3) S0 = 0.5*(1.0 - tanh(vn/u0/delta)) dongOutflow = -0.5*(ux*ux+uy*uy+uz*uz)*S0 return end c----------------------------------------------------------------------- subroutine comp_Reynolds implicit none include 'SIZE' include 'TOTAL' integer i,n,npeb real Rcyl,Rpeb real z0,z1,wd,fun_z,wt(lx1*ly1*lz1*lelv) real zz,dis,glsum,glsc2 real wght_vol,area_peb,area_cyl1,area_cyl2,vz_bar,Dh1,Dh2 npeb = 1568 Rpeb = 1.0 Rcyl = 13.858 ! min/max of pebble centers are -1.05733229e+01 1.15912238e+01 ! f ------\ ! \ ! \----- ! o ) . ! P 1 | wd | wd = 1.5 z0 = -10.573 - 1.0 - wd*0.5 z1 = 11.159 + 1.0 + wd*0.5 n=lx1*ly1*lz1*nelv do i=1,n zz=zm1(i,1,1,1) if (zz>z1) then dis = -abs(zz-z0) elseif (zz>z0) then dis = min(abs(z1-zz),abs(zz-z0)) else dis = -abs(z0-zz) endif fun_z = tanh(dis*wd) wt(i) = fun_z * bm1(i,1,1,1) enddo wght_vol = glsum(wt,n) area_peb = npeb*4.0*pi*Rpeb*Rpeb area_cyl1 = 2.0*Rcyl*pi*(z1-z0) area_cyl2 = 2.0*Rcyl*pi*(z1-z0+wd) vz_bar=glsc2(wt,vz,n)/wght_vol Dh1 = 4.0*wght_vol/(area_peb+area_cyl1) Dh2 = 4.0*wght_vol/(area_peb+area_cyl2) if(nid.eq.0)then write(*,1)z0,z1,wd,wght_vol write(*,2)area_peb,area_cyl1,area_cyl2 write(*,3)time,vz_bar,Dh1,Dh2 endif 1 format('CompRe, domain',1p4e11.3) 2 format('CompRe, s_area',1p3e11.3) 3 format('CompRe, vz/Dh ',1p4e11.3) return end c-----------------------------------------------------------------------