[Nek5000-users] Rayleigh-Benard convection (RBC) with free slip boundary conditions
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Wed Feb 12 08:05:23 CST 2014
Hi Prosenjit,
There are several 2D RB examples that would be a good
starting point. For these, there is also a companion
document here:
http://www.mcs.anl.gov/~fischer/nek5000/examples.pdf
(See section 3).
Paul
On Wed, 12 Feb 2014, nek5000-users at lists.mcs.anl.gov wrote:
> Hi Neks,
> I am new to nek5000. I am trying to use nek for studying RBC problem
> with free slip boundary conditions.
> But in the so far tried runs I got only zero values of total energy as well
> as Nusselt number which is not expected. Can anyone please help me in this
> regard.
>
> Following are the important files which I used for the test:
>
> ***********.box file***********
> input.rea
> 3
> 2
> Box
> -9 -9 3
> 0 9 1
> 0 9 1
> 0 0.2 0.8 1.0
> P ,P ,P ,P ,SYM,SYM
> P ,P ,P ,P ,SYM,SYM
>
> ***************.usr file******************
>
> subroutine rayleigh_const
>
> include 'SIZE'
> include 'INPUT'
>
> common /rayleigh_r/ rapr,ta2pr
>
> Pr = param(2)
> eps = param(75)
> Rc = param(76)
> Ta2 = param(77)
> Ra = Rc*(1.0+eps)
>
> rapr = ra*pr
> ta2pr = ta2*pr
>
> 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
>
> ffx = 0.0 ! uy*Ta2Pr
> ffy = 0.0 !- ux*Ta2Pr
> ffz = temp*rapr
> c write(6,*) ffy,temp,rapr,'ray',ieg
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine userq (ix,iy,iz,ieg)
> include 'SIZE'
> include 'TOTAL'
> include 'NEKUSE'
>
> qvol = 0.0
> 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.0 !sin(2.8284*x)*cos(3.1416*z)
> uy = 0.0 !sin(2.8284*x)*cos(3.1416*z)
> uz = 0.0 !sin(2.8284*x)*sin(3.1416*z)
> temp = 1.0 - z
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine useric (ix,iy,iz,ieg)
> include 'SIZE'
> include 'TOTAL'
> include 'NEKUSE'
> integer idum
> save idum
> data idum /99/
>
> c ran = ran1(idum)
>
> c The totally ad-hoc random number generator below is preferable
> c to the above for the simple reason that it gives the same i.c.
> c independent of the number of processors, which is important for
> c code verification.
>
> ran = 2.e4*(ieg+x*sin(y)) + 1.e4*ix*iy + 1.e5*ix*iz
> ran = 1.e3*sin(ran)
> ran = 1.e3*sin(ran)
> ran = cos(ran)
> amp = 0.1
>
> temp = 1-z + ran*amp*(1-z)*z*x*(9-x)*y*(9-y)
>
> ux=sin(2.8284*x)*cos(3.1416*z)
> uy=sin(2.8284*y)*cos(3.1416*z)
> uz=-2.8284/3.1416*cos(2.8284*x)*sin(3.1416*z)
> &-2.8284/3.1416*cos(2.8284*y)*sin(3.1416*z)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine usrdat
> return
> end
> c-----------------------------------------------------------------------
> subroutine usrdat3
> return
> end
> c-----------------------------------------------------------------------
> subroutine usrdat2
> include 'SIZE'
> include 'TOTAL'
>
> common /rayleigh_r/ rapr,ta2pr
>
> call rayleigh_const
>
>
> param(66) = 4
> param(67) = 4
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine userchk
> include 'SIZE'
> include 'TOTAL'
> parameter(lt=lx1*ly1*lz1*lelv)
>
> common /scrns/ tz(lx1*ly1*lz1*lelt)
> common /rayleigh_r/ rapr,ta2pr
> common /rayleigh_c/ Ra,Ra1,Ra2,Ra3,Prnum,Ta2,Ek1,Ek2,Ek3,ck
>
>
> real Ek0,Ek,t12,t23
> real Ey,Ex,Ez,Et,Ewt,tx,ty,tz
>
> save Ek0,Ek,t12,t23
> ra=param(76)
> prnum=param(2)
> if (nid.eq.0.and.istep.eq.0) open(unit = 79, file = 'glob.dat',
> $ status = 'new')
> n = nx1*ny1*nz1*nelv
> Ek0 = Ek
>
> Ewt = (glsc3(vz,t,bm1,n)/volvm1)
> Ey = (glsc3(vy,vy,bm1,n)/volvm1)
> Ex = (glsc3(vx,vx,bm1,n)/volvm1)
> Ez = (glsc3(vz,vz,bm1,n)/volvm1)
> Et = (glsc3(t,t,bm1,n)/volvm1)
> Ek = Ex+Ey+Ez
>
> sigma = 1.e-4
> de = abs(Ek-Ek0)/dt
>
> c if (nid.eq.0) write(79,6) istep,time,ra,prnum,Ek,Ex,Ey,Ez,
> c $ Et,Ewt
> c 6 format(i7,1p9e13.5)
>
> c n = nx1*ny1*nz1*nelv
> umax = glmax(vx,n)
> vmax = glmax(vy,n)
> wmax = glmax(vz,n)
>
> if (istep.eq.0) then
> ifxyo = .true. ! For VisIt
> do i=1,n
> tz(i) = t(i,1,1,1,1) + zm1(i,1,1,1) - 1.0
> enddo
> call outpost(vx,vy,vz,pr,tz,' ')
> else
> ifxyo = .false.
> endif
> c if (nid.eq.0) write(79,6)
> if (nid.eq.0) write(79,1)istep,time,ra,prnum,umax,vmax,wmax,
> $ Ex,Ey,Ez,Ek,Et,Ewt
> 1 format(i9,1p12e14.6)
>
> return
> end
>
> ****************.rea file ******************************
> ****** PARAMETERS *****
> 2.60000 NEKTON VERSION
> 3 DIMENSIONAL RUN
> 103 PARAMETERS FOLLOW
> 1.00000 p1 DENSITY
> 0.20000 p2 VISCOS
> 0.
> 0.
> 0.
> 0.
> 1.00000 p7 RHOCP
> 1.00000 p8 CONDUCT
> 0.
> 0. p10 FINTIME
> 75000.00 p11 NSTEPS
> -.0050000 p12 DT
> 0. p13 IOCOMM
> 0. p14 IOTIME
> 1000.000 p15 IOSTEP
> 0. p16 PSSOLVER
> 0.
> 0.250000E-01 p18 GRID
> -1.00000 p19 INTYPE
> 4.0000 p20 NORDER
> 0.000000E-06 p21 DIVERGENCE
> 0.000000E-08 p22 HELMHOLTZ
> 0 p23 NPSCAL
> 1.000000E-03 p24 TOLREL
> 1.000000E-03 p25 TOLABS
> 2.01000 p26 COURANT/NTAU
> 2.00000 p27 TORDER
> 0. p28 TORDER: mesh velocity (0: p28=p27)
> 0. p29 magnetic visc if > 0, = -1/Rm if < 0
> 0. p30 > 0 ==> properties set in uservp()
> 0. p31 NPERT: #perturbation modes
> 0. p32 #BCs in re2 file, if > 0
> 0.
> 0.
> 0.
> 0.
> 0.
> 0.
> 0.
> 0.
> 0. p41 1-->multiplicative SEMG
> 0. p42 0=gmres/1=pcg
> 0. p43 0=semg/1=schwarz
> 0. p44 0=E-based/1=A-based prec.
> 0. p45 Relaxation factor for DTFS
> 0. p46 reserved
> 0. p47 vnu: mesh matieral prop
> 0.
> 0.
> 0.
> 0.
> 0. p52 IOHIS
> 0.
> 0. p54 1,2,3-->fixed flow rate dir=x,y,z
> 0. p55 vol.flow rate (p54>0) or Ubar (p54<0)
> 0.
> 0.
> 0.
> 0. p59 !=0 --> full Jac. eval. for each el.
> 0. p60 !=0 --> init. velocity to small nonzero
> 0.
> 0. p62 >0 --> force byte_swap for output
> 0. p63 =8 --> force 8-byte output
> 0. p64 =1 --> perturbation restart
> 0. p65 #iofiles (eg, 0 or 64); <0 --> sep. dirs
> 4.00000 p66 output : <0=ascii, else binary
> 4.00000 p67 restart: <0=ascii, else binary
> 0. p68 iastep: freq for avg_all
> 0.
> 0.
> 0.
> 0.
> 0.
> 0. p74 verbose Helmholtz
> 0. p75 epsilon for RB criticality (in .usr)
> 1000. p76 Rayleigh number (in .usr)
> 0.
> 0.
> 0.
> 0.
> 0.
> 0.
> 0.
> 0. p84 !=0 --> sets initial timestep if p12>0
> 0. p85 dt ratio if p84 !=0, for timesteps>0
> 0. p86 reserved
> 0.
> 0.
> 0.
> 0.
> 0.
> 0.
> 40.0000 p93 Number of previous pressure solns saved
> 5.00000 p94 start projecting velocity after p94 step
> 5.00000 p95 start projecting pressure after p95 step
> 0.
> 0.
> 0.
> 3.00000 p99 dealiasing: <0--> off/3--> old/4--> new
> 0.
> 0. p101 No. of additional filter modes
> 1.00000 p102 Dump out divergence at each time step
> .000 p103 weight of stabilizing filter (.01)
> 4 Lines of passive scalar data follows2 CONDUCT; 2RHOCP
> 1.00000 1.00000 1.00000 1.00000 1.00000
> 1.00000 1.00000 1.00000 1.00000
> 1.00000 1.00000 1.00000 1.00000 1.00000
> 1.00000 1.00000 1.00000 1.00000
> 13 LOGICAL SWITCHES FOLLOW
> T IFFLOW
> T IFHEAT
> T IFTRAN
> T T F F F F F F F F F IFNAV & IFADVC (convection in P.S.
> fields)
> F F T T T T T T T T T T IFTMSH (IF mesh for this field is T
> mesh)
> F IFAXIS
> F IFSTRS
> F IFSPLIT
> F IFMGRID
> F IFMODEL
> F IFKEPS
> F IFMVBD
> F IFCHAR
> 8.00000 8.00000 -0.500000 -4.00000
> XFAC,YFAC,XZERO,YZERO
> *** MESH DATA ***
> 243 3 243 NEL,NDIM,NELV
>
> *****************SIZE *************************************
> C Dimension file to be included
> C
> C HCUBE array dimensions
> C
> parameter (ldim=3)
> parameter (lx1=12,ly1=lx1,lz1=lx1,lelt=31,lelv=lelt)
> parameter (lxd=18,lyd=lxd,lzd=lxd)
> parameter (lelx=1,lely=1,lelz=1)
>
> parameter (lzl=3 + 2*(ldim-3))
>
> parameter (lx2=lx1-2)
> parameter (ly2=ly1-2)
> parameter (lz2=lz1-2)
> parameter (lx3=lx1)
> parameter (ly3=ly1)
> parameter (lz3=lz1)
>
> parameter (lp = 8)
> parameter (lelg = 243)
> c
> c parameter (lpelv=lelv,lpelt=lelt,lpert=3) ! perturbation
> c parameter (lpx1=lx1,lpy1=ly1,lpz1=lz1) ! array sizes
> c parameter (lpx2=lx2,lpy2=ly2,lpz2=lz2)
> c
> parameter (lpelv=1,lpelt=1,lpert=1) ! perturbation
> parameter (lpx1=1,lpy1=1,lpz1=1) ! array sizes
> parameter (lpx2=1,lpy2=1,lpz2=1)
> c
> c parameter (lbelv=lelv,lbelt=lelt) ! MHD
> c parameter (lbx1=lx1,lby1=ly1,lbz1=lz1) ! array sizes
> c parameter (lbx2=lx2,lby2=ly2,lbz2=lz2)
> c
> parameter (lbelv=1,lbelt=1) ! MHD
> parameter (lbx1=1,lby1=1,lbz1=1) ! array sizes
> parameter (lbx2=1,lby2=1,lbz2=1)
>
> C LX1M=LX1 when there are moving meshes; =1 otherwise
> parameter (lx1m=1,ly1m=1,lz1m=1)
> parameter (ldimt= 3) ! 3 passive scalars + T
> parameter (ldimt1=ldimt+1)
> parameter (ldimt3=ldimt+3)
> c
> c Note: In the new code, LELGEC should be about sqrt(LELG)
> c
> PARAMETER (LELGEC = 1)
> PARAMETER (LXYZ2 = 1)
> PARAMETER (LXZ21 = 1)
>
> PARAMETER (LMAXV=LX1*LY1*LZ1*LELV)
> PARAMETER (LMAXT=LX1*LY1*LZ1*LELT)
> PARAMETER (LMAXP=LX2*LY2*LZ2*LELV)
> PARAMETER (LXZ=LX1*LZ1)
> PARAMETER (LORDER=3)
> PARAMETER (MAXOBJ=4,MAXMBR=LELT*6)
> PARAMETER (lhis=100) ! # of pts a proc reads from hpts.in
> ! Note: lhis*np > npoints in hpts.in
> C
> C Common Block Dimensions
> C
> PARAMETER (LCTMP0 =2*LX1*LY1*LZ1*LELT)
> PARAMETER (LCTMP1 =4*LX1*LY1*LZ1*LELT)
> C
> C The parameter LVEC controls whether an additional 42 field arrays
> C are required for Steady State Solutions. If you are not using
> C Steady State, it is recommended that LVEC=1.
> C
> PARAMETER (LVEC=1)
> C
> C Uzawa projection array dimensions
> C
> parameter (mxprev = 20)
> parameter (lgmres = 20)
> C
> C Split projection array dimensions
> C
> parameter(lmvec = 1)
> parameter(lsvec = 1)
> parameter(lstore=lmvec*lsvec)
> c
> c NONCONFORMING STUFF
> c
> parameter (maxmor = lelt)
> C
> C Array dimensions
> C
> COMMON/DIMN/NELV,NELT,NX1,NY1,NZ1,NX2,NY2,NZ2
> $,NX3,NY3,NZ3,NDIM,NFIELD,NPERT,NID
> $,NXD,NYD,NZD
>
> c automatically added by makenek
> parameter(lxo = lx1) ! max output grid size (lxo>=lx1)
>
> c automatically added by makenek
> parameter(lpart = 1 ) ! max number of particles
>
> c automatically added by makenek
> parameter (lfdm=0) ! == 1 for fast diagonalization method
>
> c automatically added by makenek
> integer ax1,ay1,az1,ax2,ay2,az2
> parameter (ax1=lx1,ay1=ly1,az1=lz1,ax2=lx2,ay2=ly2,az2=lz2) ! running
> averages
>
> c automatically added by makenek
> parameter (lxs=1,lys=lxs,lzs=(lxs-1)*(ldim-2)+1) !New Pressure
> Preconditioner
>
>
>
> With thanks and regards
> Prosenjit
>
More information about the Nek5000-users
mailing list