[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