[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 07:33:58 CST 2014
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20140212/04bee98d/attachment-0001.html>
More information about the Nek5000-users
mailing list