c----------------------------------------------------------------------- C C USER SPECIFIED ROUTINES: C C - boundary conditions C - initial conditions C - variable properties C - forcing function for fluid (f) C - forcing function for passive scalar (q) C - general purpose routine for checking errors etc. C c----------------------------------------------------------------------- subroutine uservp (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' C UDIFF =0. UTRANS=0. return end c----------------------------------------------------------------------- subroutine userf (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' EQUIVALENCE (prandtl,param(2)) common/ra_pr/rapr,Ta2_pr C FFX = 0.0 FFY = 0.0 FFZ = 0.0 bouy = (temp)*rapr if (if3d) then ffx = uy*Ta2_pr !the coriolis force terms ffy =-ux*Ta2_pr ffz = bouy elseif (ifaxis) then ffx = -bouy ffy = 0.0 ffz = 0.0 else ! 2D ffx = 0.0 ffy = bouy ffz = 0.0 endif return end c----------------------------------------------------------------------- subroutine userq (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' C QVOL = 0.0 c SOURCE = 0.0 C########## adding a reaction term to the first passive scalar c ifield equal 1 is velocity, 2 is temperature if (ifield.eq.3) then ! passive scalar 1 ps1 = ps(1) gamma = 1 qvol = gamma*ps1*(1-ps1) elseif (ifield.eq.4) then ! passive scalar 2 ps2 = ps(2) gamma = 1 qvol = gamma*ps2*(1-ps2) elseif (ifield.eq.5) then ! passive scalar 3 ps3 = ps(3) gamma = 1 qvol = gamma*ps3*(1-ps3) endif return end c----------------------------------------------------------------------- subroutine userchk include 'SIZE' include 'TOTAL' c parameter (lt=lx1*ly1*lz1*lelt) common/store/ pso(lx1,ly1,lz1,lelt),psn(lx1,ly1,lz1,lelt) common/dc_dt/ dcdt(lx1,ly1,lz1,lelt) real dt dt=param(12) if(istep.gt.0) then psn = t(:,:,:,:,2) pso = vxlag(:,:,:,:,4) dcdt = (psn-pso)/dt c print*,dcdt endif return end c----------------------------------------------------------------------- subroutine userbc (ix,iy,iz,iside,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' ux=0.0 uy=0.0 uz=0.0 c temp = 1-z + bias c arg = -istep/.005 c bias = x*0.2*exp(arg) ! use this to break axisymmetry c#### this is the usual c temp = 1-z + bias ! early in the simulation, only temp = 1-z c#### this is a hot sidewall c if (z .lt. 1.0 .and. z .gt. 0.0) then c temp = 1.0 c endif return end c----------------------------------------------------------------------- c c this makes an initial condition of random thermal perturbations c subroutine useric (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' integer idum save idum data idum /999/ c the velocity field if (ifield.eq.1) then ux = 0 uy = 0 uz = 0 elseif (ifield.eq.2) then c the temperature field temp = 1-z + 0.01*ran2(idum) c temp = 1-z + 0.1*sin(3.5*x)*sin(pi*z) c the passive scalar fields elseif (ifield.eq.3) then c initialize the passive scalars temp = exp(-(x**2+y**2)/(2*0.5**2)) elseif (ifield.eq.4) then c initialize the passive scalars c temp = sin(2*pi*(sqrt(x**2+y**2))) temp = exp(-(x**2+y**2)/(2*0.5**2)) elseif (ifield.eq.5) then c initialize the passive scalars c temp = 0.5 temp = exp(-(x**2+y**2)/(2*0.5**2)) endif c ux = 0 c uy = 0 c uz = 0 c temp = 1-z + 0.01*ran2(idum) cc temp = 1-z + 0.1*sin(3.5*x)*sin(pi*z) return end c----------------------------------------------------------------------- subroutine usrdat return end c----------------------------------------------------------------------- subroutine usrdat2 include 'SIZE' include 'TOTAL' c c This is where geometry modifications can take place c#### this was suggested by stefan -- i.e. not to hardwire the c file format but to set this in the rea file c param(66) = 4 c param(67) = 4 c try this mrp 07/14/16 param(66) = 6 param(67) = 6 return end c----------------------------------------------------------------------- subroutine usrdat3 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 c----------------------------------------------------------------------- subroutine userchk_flux include 'SIZE' include 'TOTAL' c JDS DEBUG 6/6/2008 c This routine is for post-processing data, such as computing c integrated heat fluxes, etc. c common /SCRNS/ dtdx(lx1,ly1,lz1,lelt) $ , dtdy(lx1,ly1,lz1,lelt) $ , dtdz(lx1,ly1,lz1,lelt) $ , w (lx1,ly1,lz1,lelt) common /chkcmnr/ atime,timel,flux(2),favg(2),frms(2) common /chkcmni/ iesurf(0:2*lelt,2),ifsurf(0:2*lelt,2) c integer icalld save icalld data icalld /0/ CHARACTER*3 CB1B c if (icalld.eq.0) then icalld = icalld + 1 atime = 0. timel = time call rzero(favg,2) call rzero(frms,2) c c c c c Find list of faces to know where to c compute the integrated fluxes c call find_srf(iesurf(0,1),ifsurf(0,1),5) ! bottom face call find_srf(iesurf(0,2),ifsurf(0,2),6) ! top face c endif if (istep.le.0) return c c c Compute surface fluxes c call dudxyz(dtdx,t,rxm1,sxm1,txm1,jacm1,2,1) call dudxyz(dtdy,t,rym1,sym1,tym1,jacm1,2,2) if (if3d) call dudxyz(dtdz,t,rzm1,szm1,tzm1,jacm1,2,3) c do k = 1,2 flux(k) = 0.0 do isurf = 1,iesurf(0,k) ie = iesurf(isurf,k) iface = ifsurf(isurf,k) CB1B = CBC(IFACE,IE,1) if (CB1B.EQ.'E ') then c do nothing else call surface_flux(dq,dtdx,dtdy,dtdz,ie,iface,w) flux(k) = flux(k) + dq end if c write(6,*) isurf,' dq',ie,iface,k,dq enddo enddo c Sum over all processors call gop(flux,w,'+ ',2) c c c c Compute time average, rms c dtime = time - timel atime = atime + dtime if (nid.eq.0.and.atime.ne.0..and.dtime.ne.0.) then beta = dtime/atime alpha = 1.-beta do k=1,2 favg(k) = alpha*favg(k) + beta*flux(k) frms(k) = alpha*frms(k) + beta*flux(k)*flux(k) c frrms = frms(k) - favg(k)**2 frrms = sqrt(frrms) c write(6,1) istep,k,time,atime,flux(k),favg(k),frrms enddo elseif (nid.eq.0) then do k=1,2 write(6,1) istep,k,time,atime,flux(k),favg(k),frms(k) enddo endif 1 format(i6,i2,' flx',1p5e14.6) timel = time c return end c----------------------------------------------------------------------- c----------------------------------------------------------------------- subroutine find_srf(iesurf,ifsurf,inface) c c Find list of surfaces over which the flux should be computed c c Number of such surfaces is returned in iesurf(0). c c include 'SIZE' include 'TOTAL' c integer iesurf(0:lelt),ifsurf(0:lelt) c nsurf = 0 nfaces = 2*ndim do ie=1,nelv nsurf = nsurf+1 iesurf(nsurf) = ie ifsurf(nsurf) = inface enddo iesurf(0) = nsurf ifsurf(0) = nsurf c write(6,*) nid,' num surfaces',nsurf return end c----------------------------------------------------------------------- c----------------------------------------------------------------------- function ran2(delta) c real*8 F7,T PARAMETER ( F7=78125., T30=1073741824.) c save t c integer icalld save icalld data icalld/0/ C C INITIALIZATION: C Generate PSEUDO-RANDOM (0., 1.) DATA C USING A RANDOM NUMBER GENERATOR BASED ON THE RECURSION C X(N+1) = 5**7 * X(N) (MOD 2**30) C THIS RECURSION WILL GENERATE 2**28 (APPROX. 268 MILLION) NUMBERS C BEFORE REPEATING. FOR THIS SCHEME TO WORK PROPERLY, THE HARDWARE C MULTIPLY OPERATION MUST BE CORRECT TO 47 BITS OF PRECISION. C if (icalld.eq.0) T = F7 / T30 icalld = 1 T = MOD (F7 * T, 1.) ran2 = delta*t return end c-----------------------------------------------------------------------