c----------------------------------------------------------------------- c c A short example of Rayleigh-Benard convection: c c Parameters are set in routine rayleigh_const for convenience. c c With this nondimensionalization, set rho==1 (parameter p1 in .rea c file) and visocity (p2) to be the desired Prandtl number. c c Rayleigh number is set as Ra = Rc*(1+eps), Rc=p76, eps=p75. c c The buoyancy is ffy = Ra Pr T, where T is determined by c boundary and initial conditions. c c Critical Rayleigh number is around 1707.762 c (Somehow I was recalling 1734, but that appears to be for a c particular geometric configuration considered by Laurette Tuckerman c & Dwight Barkley) c c GEOMETRY: c c There are two primary cases, ray1.box and ray2.box. c The former specifies 10 elements in x, the latter only 9, c both for a 9x1 domain. c c NOTES: c c A time trace of (u,v)_max vs t is output to the logfile. See userchk. c c Be careful about selecting an even number of elements in x c as it appears that the RB system likes to lock onto the grid spacing c and give a number of rolls that matches the number of elements, if the c elements have order-unity aspect ratio, as in the present case. c Thus, in the case, the 9 element mesh is likely to be more faithful c to the linear stability theory, at least for modest polynomial orders c of lx1=12. c c It appears that one cannot realize Courant conditions of CFL ~ 0.5 c with these cases because of the explicit Boussinesq treatment. c The given value dt=.02 is stable with lx1=12. c c c----------------------------------------------------------------------- subroutine rayleigh_const include 'SIZE' include 'INPUT' common /rayleigh_r/ rapr,ta2pr common /myinputs/ ra,pra,rai,tramp,tbegav,extol,dnsave,nout,dnout, $ horout,tout,ostep,dtout common /myconst/ q0,scalet integer dnsave,nout,horout,dnout,ostep Pra = param(2) eps = param(75) c Rc = param(76) c Ra = Rc*(1.+eps) Ra = param(76) Ta2 = param(77) Rai = param(121) tramp = param(122) tbegav = param(123) extol = param(124) dnsave = param(125) nout = param(126) dnout = param(127) horout = param(128) tout = param(129) ostep = param(130) dtout = param(131) scalet = param(132) rapr = rai*pra ta2pr = ta2*pra return end c----------------------------------------------------------------------- subroutine uservp (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' udiff = 0 utrans = 0 return end c----------------------------------------------------------------------- subroutine userf (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' common /rayleigh_r/ rapr,ta2pr buoy = temp*rapr if (if3d) then ffx = uy*Ta2Pr ffy = - ux*Ta2Pr ffz = buoy else if (ifaxis) then ffx = -buoy ffy = 0. else ffx = 0. ffy = buoy end if c write(6,*) ffy,temp,rapr,'ray',ieg return end c----------------------------------------------------------------------- subroutine userq (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' common /myconst/ q0,scalet q0 = 1 if (ifield.eq.2) then qvol = (temp+2*(1-y))**4*ps(1)**2 else if (ifield.eq.3) then qvol = -(temp+2*(1-y))**4*ps(1)**2 end if source = 0.0 return end c----------------------------------------------------------------------- subroutine userbc (ix,iy,iz,iside,ieg) include 'SIZE' include 'TSTEP' include 'INPUT' include 'NEKUSE' common /rayleigh_r/ rapr,ta2pr ux=0. uy=0. uz=0. if (ifield.eq.2) then temp = 0. flux = 0. else if (ifield.eq.3) then temp = 0. flux = 0. end if c if (y.gt.0.5) then c flux = -1 c else c flux = 1 c end if c if (if3d) then c temp = 1-z c else if (ifaxis) then ! domain is on interval x in [-1,0] c temp = 1.+x c else ! 2D: domain is on interval y in [0,1] c temp = 1.-y c end if return end c----------------------------------------------------------------------- subroutine useric (ix,iy,iz,ieg) include 'SIZE' include 'TOTAL' include 'NEKUSE' integer idum save idum data idum /99/ c ran = 2.e7*(ieg+x*sin(y)) + 1.e6*ix*iy + 1.e7*ix c ran = 1.e9*sin(ran) c ran = 1.e9*sin(ran) c ran = cos(ran) c ran = ran1(idum) c amp = 1E-6 if (ifield.eq.2) then temp = 0 else if (ifield.eq.3) then temp = 2 end if c temp = 1-y + ran*amp*(1-y)*y*x*(2-x) c temp = .5*y*(1-y) + ran*amp*(1-y)*y ux=0.0 uy=0.0 uz=0.0 return end c----------------------------------------------------------------------- subroutine usrdat return end c----------------------------------------------------------------------- subroutine usrdat3 return end c----------------------------------------------------------------------- subroutine usrdat2 include 'SIZE' include 'TOTAL' common /rayleigh_r/ rapr,ta2pr common /myinputs/ ra,pra,rai,tramp,tbegav,extol,dnsave,nout,dnout, $ horout,tout,ostep,dtout call rayleigh_const param(66) = 4 param(67) = 4 return end c----------------------------------------------------------------------- subroutine userchk include 'SIZE' include 'TOTAL' include 'ZPER' include 'AVG' common /scrns/ tz(lx1*ly1*lz1*lelt) common /achk/ yh(ly1*lely) $ , w1(ly1*lely) $ , w2(ly1*lely) c $ , vha(ly1*lely) $ , th(ly1*lely) $ , ch(ly1*lely) c $ , tha(ly1*lely) common /rayleigh_r/ rapr,ta2pr common /myinputs/ ra,pra,rai,tramp,tbegav,extol,dnsave,nout,dnout, $ horout,tout,ostep,dtout real Tav,Cav,J,Ek,at_av,aC_av,aJ,aEk,cumJ,cumEk,cumt,cumC real Ektol,tavtol,Jtol real oldEk(0:400),oldt(0:400),oldJ(0:400) integer icount,isave,avgstp,endstp,dnsave,dnout,nout,horout,ostep logical avbeg save oldEk,oldt,oldJ,oldC save avbeg data avbeg /.false./ save icount data icount /0/ save isave data isave /0/ save avgstp data avgstp /1E8/ save endstp data endstp /1E8/ nelx=lelx nely=lely nelz=lelz ny1=ly1 ny = ny1*nely c write(6,*) 'nid= ',nid c If ramp-up over, give Ra its final value if (time.ge.tramp) rapr=ra*pra c Calculate horizontal average of time average call avg_all call planar_average_y(yh,ym1,w1,w2) if (mod(istep,horout).eq.0) then icount=icount+1 write(6,*) istep,time,' tpts2=',t(2,2,1,2,1),t(2,2,1,2,2) call planar_average_y_ps(th,t,w1,w2,1) call planar_average_y_ps(ch,t,w1,w2,2) c write(6,*) istep,time,' th=',th,' ch=',ch c call planar_average_y(tha,tavg,w1,w2) c write(6,2) istep,time,' at_av= ',at_av,' call out1d_2(yh,th,ch,'yTC1d',ny,icount) c call out1d(yh,ch,'ych1d',ny,icount) c n = nx1*ny1*nz1*nelt c at_int = glsc2(tavg,bm1,n) c at_av = at_int/voltm1 c aEk = .5*(glsc3(uavg,bm1,uavg,n)+glsc3(vavg,bm1,vavg,n)+ c $ glsc3(wavg,bm1,wavg,n))/voltm1 c aJ = glsc3(vavg,bm1,tavg,n)/voltm1 c if (nid.eq.0) then c write(6,2) istep,time,' at_av= ',at_av,' aEk= ',aEk, c $ ' aJ= ',aJ c end if c 2 format(i8,e12.4,a,e10.4,a,e10.4,a,e10.4) end if c Calculate maximum velocities n = nx1*ny1*nz1*nelv umax = glmax(vx,n) vmax = glmax(vy,n) c Output initial fields and maximum velocities ifxyo = .true. ! For VisIt if (istep.gt.iostep) ifxyo = .false. if (istep.le.0) then do i=1,n tz(i) = t(i,1,1,1,1)+ym1(i,1,1,1) - 1. enddo call outpost2(vx,vy,vz,pr,t,2,' ') end if if (nid.eq.0) write(6,1) istep,time,umax,vmax 1 format(i9,1p3e14.6,' umax') c Calculate and write volume averages n = nx1*ny1*nz1*nelt Tav = glsc2(t,bm1,n)/voltm1 n = nx1*ny1*nz1*nelfld(2) Cav = glsc2(bm1,t(1,1,1,1,2),n)/volfld(2) Ek = .5*(glsc3(vx,bm1,vx,n)+glsc3(vy,bm1,vy,n)+ $ glsc3(vz,bm1,vz,n))/voltm1 J = glsc3(vy,bm1,t,n)/voltm1 if (nid.eq.0) write(6,3) istep,time,' Tav= ',Tav,' Cav= ',Cav, $ ' Ek= ',Ek,' J= ',J 3 format(i8,e12.4,a,e10.4,a,e10.4,a,e10.4,a,e10.4) c Set step at which cumulative averaging begins if (time.ge.tbegav.and.(.not.avbeg)) then avbeg=.true. avgstp=istep c if (nid.eq.0) write(6,*) 'avgstp= ',avgstp write(6,*) 'avgstp= ',avgstp,nid end if c if (nid.eq.0) then c write(*,*) 'tbegav=',tbegav,avbeg c write(*,*) 'avgstp=',avgstp,istep c end if c Cumulative averages if (time.ge.tbegav) then call runtimeavg(cumJ,J,1,avgstp,1000,'cumJ ') call runtimeavg(cumt,Tav,2,avgstp,1000,'cumt ') call runtimeavg(cumEk,Ek,3,avgstp,1000,'cumEk') end if c Store cumulative averages and check tolerances if (time.ge.tbegav.and.mod(istep,dnsave).eq.0) then oldEk(isave) = cumEk oldt(isave) = cumt oldJ(isave) = cumJ c if (nid.eq.0) write(*,*) 'isave= ',isave write(6,*) 'isave= ',isave,nid if (isave.ge.2.and.mod(isave,2).eq.0) then c if (nid.eq.0) write(*,*) 'oldEk(i/2)= ',oldEk(isave/2) Ektol = abs((cumEk-oldEk(isave/2))/cumEk) tavtol = abs((cumt-oldt(isave/2))/cumt) Jtol = abs((cumJ-oldJ(isave/2))/cumJ) if (nid.eq.0) write(*,*) istep,' tols: ',Ektol,tavtol,Jtol if (nid.eq.0) write(*,*) istep,'maxtol= ', $ max(Ektol,tavtol,Jtol) c if (max(Ektol,tavtol,Jtol).lt.extol) then cc if (nid.eq.0) write(*,*) 'tolerances met' c write(*,*) istep,'tolerances met',nid c if (endstp.ge.1E8) c $ endstp = istep - mod(istep,dnout) + dnout*nout c write(*,*) 'endstp= ',endstp c end if end if isave = isave+1 end if c Output end frames c Exit if end frames done c if (istep.ge.endstp) then c write (6,*) 'frames done, calling exitt' c call gsync c call exitt c end if c Output end frames if (time.ge.tout.and.time.le.(tout+dtout).and. $ mod(istep,ostep).eq.0) then call outpost(vx,vy,vz,pr,t,' ') end if return end c----------------------------------------------------------------------- subroutine runtimeavg(ay,y,j,istep1,ipostep,s5) c c Computes, stores and for ipostep!=0 prints runtime averages c of j-quantity y (along w/ y itself unless ipostep<0) c with j + 'rtavg_' + (unique) s5 every ipostep for istep>=istep1 c include 'SIZE' ! for nid & TSTEP include 'TSTEP' ! for istep & time (or TOTAL) parameter (lymax=99) ! max of averages to store character*5 s5 ! unique string to append to 'rtavg_' character*14 s14 ! i3 + 'rtavg_' + s5 real a (lymax) ! run time average character*5 as5(lymax) ! j's strings s5 save a, as5 data a /lymax*0./, as5 /lymax*' '/ if (nid.eq.0.and.istep.ge.istep1) then ! otherwsie skip istep c--- if (istep1.lt.0) then istep1 = 0 write(6,*) 'Warning: istep1<0 in runtimeavg -- resetting to 0' end if if (as5(j).ne.s5.and.as5(j).ne.' ') then ! check for unique j & s5 write(6,10) j,s5,as5(j) call exitt end if iistep = istep - istep1 + 1 if (iistep.gt.1) then ! a_ii = (1-1/ii)*a_ii-1 + 1/ii*y_ii wy = 1./iistep wa = 1. - wy ay = wy*y + wa*a(j) ! runtime average a_ii else if (iistep.eq.1) then ! skip istep and < u^2 > real u(nx1,ny1,nz1,nelv) ! Input real w1(ny1,nely),w2(ny1,nely) ! work arrays integer e,eg,ex,ey,ez ny = ny1*nely call rzero(ua,ny) c call rzero(u2,ny) call rzero(w1,ny) do e=1,nelt eg = lglel(e) call get_exyz(ex,ey,ez,eg,nelx,nely,nelz) do k=1,nz1 do i=1,nx1 do j=1,ny1 yy = (1.-zgm1(j,2))/2. ! = 1 for j=1, = 0 for j=ny1 aa = yy*area(i,1,1,e) + (1-yy)*area(i,1,3,e) ! wgtd jacobian, fc 1&3 w1(j,ey) = w1(j,ey) + aa ua(j,ey) = ua(j,ey) + aa*u(i,j,k,e) c write(*,*) 'UA',ua(j,ey) c u2(j,ey) = u2(j,ey) + aa*u(i,j,k,e)**2 enddo enddo enddo enddo call gop(ua,w2,'+ ',ny) c call gop(u2,w2,'+ ',ny) call gop(w1,w2,'+ ',ny) do i=1,ny ua(i,1) = ua(i,1) / w1(i,1) ! Normalize c u2(i,1) = u2(i,1) / w1(i,1) ! Normalize enddo return end c----------------------------------------------------------------------- subroutine planar_average_y_ps(ua,u,w1,w2,tfld) c c Compute r-s planar average of quantity u() and u^2() c include 'SIZE' include 'GEOM' include 'PARALLEL' include 'WZ' include 'ZPER' include 'TSTEP' integer e,eg,ex,ey,ez,tfld real ua(ny1,nely) ! ,u2(ny1,nely) ! Output: < u > and < u^2 > real u(lx1,ly1,lz1,lelt,ldimt) ! Input real w1(ny1,nely),w2(ny1,nely) ! work arrays c write(6,*) ' lx1=',lx1,' lz1=',lz1,' lelt=',lelt,' ldimt=',ldimt c write(6,*) ' nelfld=',nelfld(tfld+1),' tfld=',tfld write(6,*) istep,time,' tpts=',u(2,2,1,2,1),u(2,2,1,2,2) ny = ny1*nely call rzero(ua,ny) c call rzero(u2,ny) call rzero(w1,ny) do e=1,nelfld(tfld+1) eg = lglel(e) call get_exyz(ex,ey,ez,eg,nelx,nely,nelz) do k=1,nz1 do i=1,nx1 do j=1,ny1 yy = (1.-zgm1(j,2))/2. ! = 1 for j=1, = 0 for j=ny1 aa = yy*area(i,1,1,e) + (1-yy)*area(i,1,3,e) ! wgtd jacobian, fc 1&3 w1(j,ey) = w1(j,ey) + aa ua(j,ey) = ua(j,ey) + aa*u(i,j,k,e,tfld) c write(*,*) 'UA',ua(j,ey) c u2(j,ey) = u2(j,ey) + aa*u(i,j,k,e,tfld)**2 enddo enddo enddo enddo call gop(ua,w2,'+ ',ny) c call gop(u2,w2,'+ ',ny) call gop(w1,w2,'+ ',ny) do i=1,ny ua(i,1) = ua(i,1) / w1(i,1) ! Normalize c u2(i,1) = u2(i,1) / w1(i,1) ! Normalize enddo 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