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