[Nek5000-users] Traveling gravity waves
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Thu Dec 6 06:38:42 CST 2012
Hi Neks,
I am trying to set up a traveling gravity wave and started from the
std.wv example.
- In the useric I prescribed the the Initial velocity by using linear
wave theory (Airy)
http://en.wikipedia.org/wiki/Airy_wave_theory#Solution_for_a_progressive_monochromatic_wave
wave amplitude
- the initial mesh was generate with genbox with periodic sides, bottom
wall and a free surface
- the surface tension was excluded by setting sigma & We = 0
- reverted density & viscosity to dimensional variables (SI Units) in
the .rea to get correct velocities
- set the gravitational force in userf ffy = -9.81*param(1)
But instead of a traveling wave the system reverts to a sort of standing
wave immediately
Do I have set the initial velocity of the moving mesh as well? And if so
are wx,wy the correct variables to do so?
Did I apply the dimensions correctly?
I attached the .rea & .usr, any help with my Problem would be greatly
appreciated.
Thanks
Arne
-------------- next part --------------
****** PARAMETERS *****
2.610000 NEKTON VERSION
2 DIMENSIONAL RUN --- Based on "leeho1a.rea"
120 PARAMETERS FOLLOW--- See Lee W. Ho 1989 MIT Ph.D. Thesis
998.0 p1 DENSITY
1 p2 VISCOS -3.e4 0 0.01
0.
0.
0.
0.
1.00000 p7 RHOCP
0.10000E-01 p8 CONDUCT
0.
0. p10 FINTIME
1000 p11 NSTEPS
0.000250 p12 DT
0. p13 IOCOMM
0. p14 IOTIME
5. p15 IOSTEP
0. p16 PSSOLVER
0.
0.50000E-01 p18 GRID
-1.00000 p19 INTYPE
5.0000 p20 NORDER
1.000000E-10 p21 DIVERGENCE
1.000000E-12 p22 HELMHOLTZ
0 p23 NPSCAL
0.100000E-01 p24 TOLREL
0.100000E-01 p25 TOLABS
0.40000 p26 COURANT/NTAU
3.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.0950000E+00 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.
1.00000 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
1.00000 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.01 p75 perturbation amplitude (in .usr)
0.
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.
0. p93 Number of previous pressure solns saved
0. p94 start projecting velocity after p94 step
0. p95 start projecting pressure after p95 step
0.
0.
0.
3.00000 p99 dealiasing: <0--> off/3--> old/4--> new
0.
1.00000 p101 No. of additional filter modes
1.00000 p102 Dump out divergence at each time step
-0.01 p103 weight of stabilizing filter (.01)
0.
0.
0.
0. p107 !=0--> add to h2 array in hlmhotz eqn
0.
0.
0.
0.
0.
0.
0.
0.
0 6. p116 !=0: x elements for fast tensor product
0 10. p117 !=0: y elements for fast tensor product
0 1. p118 !=0: z elements for fast tensor product
0.
0.
4 Lines of passive scalar data follows2 CONDUCT; 2RHOCP
-10.0000 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
F 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
T IFSTRS
F IFSPLIT
F IFMGRID
F IFMODEL
F IFKEPS
T IFMVBD
F IFCHAR
8.00000 8.00000 -0.500000 -4.00000 XFAC,YFAC,XZERO,YZERO
**MESH DATA** 1st line is X of corner 1,2,3,4. 2nd line is Y.
-120 2 120 NEL,NDIM,NELV
0 PRESOLVE/RESTART OPTIONS *****
7 INITIAL CONDITIONS *****
C Default
C Default
C Default
C Default
C Default
C Default
C Default
***** DRIVE FORCE DATA ***** BODY FORCE, FLOW, Q
4 Lines of Drive force data follow
C
C
C
C
***** Variable Property Data ***** Overrrides Parameter data.
1 Lines follow.
0 PACKETS OF DATA FOLLOW
***** HISTORY AND INTEGRAL DATA *****
0 POINTS. Hcode, I,J,H,IEL
***** OUTPUT FIELD SPECIFICATION *****
6 SPECIFICATIONS FOLLOW
T COORDINATES
T VELOCITY
T PRESSURE
F TEMPERATURE
F TEMPERATURE GRADIENT
0 PASSIVE SCALARS
***** OBJECT SPECIFICATION *****
0 Surface Objects
0 Volume Objects
0 Edge Objects
0 Point Objects
-------------- next part --------------
c-----------------------------------------------------------------------
c
c Setup for 2D free-surface standing wave problem
c
c Re = 998.0/1
c
c We = 0
c
c Fr = 1/9.81
c
c
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 /fsparam/ amp,alpha,Re,Pg,We,gamma_r,gamma_i,w_k
ffx = 0.0
c ffy = -(Re*Pg)**(-2)
c ffy = -1
! ffy = -9.81
ffy = -9.81*param(1)
ffz = 0.0
return
end
c-----------------------------------------------------------------------
subroutine userq (ix,iy,iz,ieg)
include 'SIZE'
include 'TOTAL'
include 'NEKUSE'
C
qvol = 0.0
source = 0.0
return
end
c-----------------------------------------------------------------------
subroutine userchk
include 'SIZE'
include 'TOTAL'
common /outputs/ vol0,ymx0,vy20,vy2l,total_ke0
common /fsarray/ ub(lx1,ly1,lz1,lelt),vb(lx1,ly1,lz1,lelt)
common /fsparam/ amp,alpha,Re,Pg,We,gamma_r,gamma_i,w_k
real dtmax
save dtmax
integer xxx
n = nx1*ny1*nz1*nelv
vol = glsum(bm1,n)
vy2 = glsc3(vy,bm1,vy,n)/vol
if (vy2.gt.0) vy2 = sqrt(vy2)
!Compute maximal surface elevation
ymx = glmax(ym1,n)
c !Compute total kinetic energy
xke = glsc3(vx, vx, bm1, n)
yke = glsc3(vy, vy, bm1, n)
total_ke = 0.5*(xke + yke)
!compute maximal mesh velocity
w_wx = glmax(wx,n)
!compute maximal velocity
w_vx = glmax(vx,n)
c
c Target error values for ge are asympotically ~ 1.e-4
c
c !Only first step
if (istep.eq.0) then
vol0 = vol
ymx0 = ymx
vy20 = vy2
c !Initial kinetic energy
total_ke0 = total_ke
c !Output Initial Condition
call outpost(vx,vy,vz,pr,t,'IC_')
!Only on Node 0
if (nid.eq.0) then
c Create Output file for Volume
OPEN (22,FILE='volume.txt',STATUS='REPLACE',
& ACCESS='sequential',position='append')
WRITE(22,*) 'Conservation of Volume'
CLOSE(20)
endif
else
c !Volume
volr = vol/vol0
ymxr = ymx/ymx0
vy2l = vy2/vy2l
!Only on Node 0
if (nid.eq.0) then
write(6,1) istep,time,ymx,vy2,ymxr
1 format(i5,1p4e14.6,' amp')
c Open Output file for Volume
OPEN (22,FILE='volume.txt',STATUS='OLD',ACCESS='sequential',
& position='append')
WRITE(22,*) time,'step', istep, 'delta Volume',(vol-vol0)/vol0
WRITE(22,*) 'kinetic energy' , total_ke,
& 'delta', (total_ke0-total_ke)
WRITE(22,*) 'maximal surface elevation' ,ymx,'delta',ymx-ymx0
WRITE(22,*) 'maximal x-mesh velocity' ,w_wx ,
& 'maximal x-velocity' ,w_vx
CLOSE(22)
endif
endif
vy2l = vy2
c if (istep.eq.0) dtmax = abs(param(12)) ! small steps at start
c if (istep.eq.0) dt = dtmax/100.
c dt = dt*1.01
c dt = min(dtmax,dt)
c param(12) = -dt
param(93) = 0 ! no projection for free-surface
param(94) = 0 ! no velocity projection for free-surface
param(95) = 0 ! no pressure projection for free-surface
return
end
c-----------------------------------------------------------------------
subroutine userbc (ix,iy,iz,iside,ieg)
include 'SIZE'
include 'TOTAL'
include 'NEKUSE'
common /fsparam/ amp,alpha,Re,Pg,We,gamma_r,gamma_i,w_k
! ux=0.0
! uy=0.0
! uz=0.0
! temp=0.0
c sigma = We ! surface tension = Weber number
c sigma = 1.e-8
sigma = 0
return
end
c-----------------------------------------------------------------------
subroutine useric (ix,iy,iz,eg)
include 'SIZE'
include 'TOTAL'
include 'NEKUSE'
common /fsarray/ ub(lx1,ly1,lz1,lelt),vb(lx1,ly1,lz1,lelt)
common /fsparam/ amp,alpha,Re,Pg,We,gamma_r,gamma_i,w_k
integer e,eg
c Prescribe Velocity
c ux=0.0
c uy=0.0
c 2D Simualation z = 0
c uz=0.0
! temp=0
c amp = 0.05 ! Perturbation amplitude
c amp = param(75)
c Wave_k ... etc
c w_k = 2.0
w_g = 9.81
c w_g = 1 !skalierte Gravitation
w_w = sqrt(w_g * w_k)
one = 1.
pi = 4.*atan(one)
x0 = 0.
x1 = 2.
n = nx1*ny1*nz1*nelv
xmin = x0
xmax = x1
do i=1,n ! perurbation displacement
x = xm1(i,1,1,1)
y = ym1(i,1,1,1)
! ax = w_k * 2.0 *pi*(x-xmin)/(xmax-xmin)
ax = w_k * x
disp1 = sin(ax)
disp2 = cos(ax)
! w_koeff = w_g * amp * w_k / w_w
w_koeff = amp * w_w
c yscl = 1+y ! zero when y=-1, 1 when y = 0
c vx(i,1,1,1) = amp * w_w * exp(w_k*y)*disp1
c vy(i,1,1,1) = -amp * w_w * exp(w_k*y)*disp2
c !set velocity
vx(i,1,1,1) = w_koeff * exp(w_k*y)*disp1
vy(i,1,1,1) = w_koeff * exp(w_k*y)*disp2
c !set mesh velocity
! wx(i,1,1,1) = w_koeff * exp(w_k*y)*disp1
! wy(i,1,1,1) = w_koeff * exp(w_k*y)*disp2
enddo
return
end
c-----------------------------------------------------------------------
subroutine usrdat
c Verbose timing
! iverbose = 1 ! fully verbose
! call platform_timer(iverbose)
! call exitti('timer done$',iverbose)
return
end
c-----------------------------------------------------------------------
subroutine usrdat3
return
end
c-----------------------------------------------------------------------
subroutine usrdat2
include 'SIZE'
include 'TOTAL'
include 'ZPER' ! nelx,nely for quickmv
include 'CTIMER' ! icalld
common /fsparam/ amp,alpha,Re,Pg,We,gamma_r,gamma_i,w_k
common /xyorig/ xo(lx1*ly1*lz1*lelt),yo(lx1*ly1*lz1*lelt)
character*3 cb
integer e,eg,ex,ey,ez
nelx = 6
nely = 6
param(66) = 4 ! std nek5000 output
param(67) = 4
c Re = 3.e4 ! hard-code Reynolds number
c visc = 1/Re
c param(2) = visc
c cpfld(1,1) = visc
c cpfld(1,2) = 1.0 ! density
! ! Pg = 1.1e-4 ! Gravitational Peclet number
We = 0.0
amp = param(75) ! Perturbation amplitude
c amp = 0.05 ! Perturbation amplitude
one = 1.
pi = 4.*atan(one)
x0 = 0.
x1 = 2.0 * pi
call rescale_x (xm1,x0,x1) ! x on [0,1] later *2*pi
y0 = -3.5
y1 = 0.
call rescale_x (ym1,y0,y1) ! y on [-1,0]
w_k = 1.0
n = nx1*ny1*nz1*nelv
xmin = glmin(xm1,n)
xmax = glmax(xm1,n)
icalld = icalld+1
do i=1,n ! perurbation displacement
x = xm1(i,1,1,1)
y = ym1(i,1,1,1)
! ax = w_k * 2*pi*(x-xmin)/(xmax-xmin)
ax = w_k * x
disp = sin(ax)
c yscl = 1+y ! zero when y=-1, 1 when y = 0
c ym1(i,1,1,1) = y + amp*yscl*disp
c ym1(i,1,1,1) = y + (y - y0) * amp * disp !same function
ym1(i,1,1,1) = y + (y - y0)/(-y0) * amp * disp !same function
enddo
param(59) = 1 ! force std operator evaluation for deformed geometry
c Set surface tension on free surface
! sigma = We ! must be set in userbc for cb = 'ms '
sigma = 0.0
itop = 3
if (if3d) itop = 5
do e=1,nelv
cb = cbc(itop,e,1) ! bc on top surface
if (cb.eq.'MS ') bc(4,itop,e,1) = sigma
enddo
return
end
c-----------------------------------------------------------------------
More information about the Nek5000-users
mailing list