Dear all,
I'm a beginner in nek5000, and I want to use the planar_average_r() function to achieve averaging along the z-direction(I'm using an older version). Can you please guide me on how to do this? Currently, I've obtained some values as output, but the total count of these values doesn't seem to match the total number of two-dimensional plane nodes. I have assigned nelx, nely, and nelz as the numbers of elements in the x-direction, y-direction, and z-direction, respectively (N-1, where N is the number of nodes). My goal is to obtain a three-dimensional flow field and calculate the average velocity along the z-direction (spanwise). Could you please advise me on how to modify the approach? I would greatly appreciate it if you could provide me with some suggestions!
Below is a portion of my .usr code :
c-----------------------------------------------------------------------
subroutine userchk
INCLUDE 'SIZE'
INCLUDE 'TOTAL'
INCLUDE 'ZPER'
c INCLUDE 'RESTART'
c Output vorticity
parameter(lt=lx1*ly1*lz1*lelv)
common /myexact/ud(lt),vd(lt)
common /myenrgy/ vx2(lx1,ly1,lz1,lelv),
& vy2(lx1,ly1,lz1,lelv),
& vz2(lx1,ly1,lz1,lelv)
common /myplane/ uavg(lx1*lelx),
& vavg(lx1*lelx),
& wavg(lx1*lelx),
& w1(lx1*lelx),
& w2(lx1*lelx),
& xx(lx1*lelx),
& yy(lx1*lelx)
character*80 filename(9999)
nelx = 144 ! Number of elements in x,y, and z directions.
nely = 19 ! NOTE, this may vary from one mesh to the next.
nelz = 9 !
write(6,*) ' Num :=', lelx, lely, lelz
ntot = nx1*ny1*nz1*nelv
if (nid .eq. 0) then
open(unit=199,file='file.list',form='formatted',status='old')
read(199,*) nfiles
read(199,'(A80)') (filename(i),i=1,nfiles)
close(199)
endif
call bcast(nfiles,isize)
call bcast(filename,nfiles*80)
do j = 1, nfiles
call load_fld(filename(j))
enddo
do i=1,ntot
vx2(i,1,1,1) = vx(i,1,1,1)**2
vy2(i,1,1,1) = vy(i,1,1,1)**2
vz2(i,1,1,1) = vz(i,1,1,1)**2
enddo
call planar_average_r(uavg,vx,w1,w2)
call planar_average_r(vavg,vy,w1,w2)
call planar_average_r(wavg,vz,w1,w2)
if(nid.eq.0) then
open(unit=56,file='avgeng.dat')
m = nx1*nelx
do i=1,m
write(56,3) uavg(i),vavg(i)
3 format(1p115e17.9)
enddo
close(56)
endif
call exitt
return
end
c-----------------------------------------------------------------------
subroutine planar_average_r(ua,u,w1,w2)
c
c Compute s-t planar average of quantity u()
c
include 'SIZE'
include 'GEOM'
include 'PARALLEL'
include 'WZ'
include 'ZPER'
c
real ua(nx1,nelx),u(nx1,ny1,nx1,nelv),w1(nx1,nelx),w2(nx1,nelx)
integer e,eg,ex,ey,ez
c
nx = nx1*nelx
call rzero(ua,nx)
call rzero(w1,nx)
c
do e=1,nelt
c
eg = lglel(e)
call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
c
do i=1,nx1
do k=1,nz1
do j=1,ny1
zz = (1.-zgm1(i,1))/2. ! = 1 for i=1, = 0 for k=nx1
aa = zz*area(j,k,4,e) + (1-zz)*area(j,k,2,e) ! wgtd jacobian
w1(i,ex) = w1(i,ex) + aa
ua(i,ex) = ua(i,ex) + aa*u(i,j,k,e)
enddo
enddo
enddo
enddo
c
call gop(ua,w2,'+ ',nx)
call gop(w1,w2,'+ ',nx)
c
do i=1,nx
ua(i,1) = ua(i,1) / w1(i,1) ! Normalize
enddo
return
end
c-----------------------------------------------------------------------
Best,
Dewu