[Nek5000-users] Help with planar averaging for dyn smag
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Fri Feb 18 21:07:55 CST 2011
Hi Harish,
I've added a function that will homogenize any function on mesh
one (lx1,ly1,lz1) in the z direction, assuming you know how many
elements you have in the x-y plane, as would be the case if your
mesh is generated using n2to3, say.
To use it, you must set lelx and lely in SIZE such that the
product lelx*lely >= nelxy. And you similarly need to define
nelx and nely values. Here are some sample lines from my userchk
routine:
call outpost(vx,vy,vz,pr,t,' ') ! Output the orginal field
nelx = 8 ! Pick an nelx,nely pair
nely = 8 ! such that nelx*nely matches my domain
call z_distribute(vx) ! Homogenize vx in the z-direction
call z_distribute(vy) ! Homogenize vy in the z-direction
call z_distribute(vz) ! Homogenize vz in the z-direction
call z_distribute(pr) ! Homogenize pr (only for Pn-Pn!!)
call outpost(vx,vy,vz,pr,t,' ') ! output homogenized field and quit
call exitti('quit after z-distribute$',nelx)
I've committed "z_distribute" to the svn repo.
Paul
On Mon, 14 Feb 2011, Harish Kanchi wrote:
> Hello Paul,
>
> I have been working on setting up a case for dyn smag. It is the same case
> that ran with Pn/Pn and smaller dt.
>
> It has NEL = 656, made up for 4 blocks. I build the 3D mesh with genbox.
> Here is the genbox file used to produce the mesh.
>
> In my case the homogeneous direction is the z-direction. Here is the
> subroutine I am using for planar averaging. But it does not work!!
>
> Any help would be appreciated.
>
> Thanks,
>
> Harish.
>
> *************************************************************************
> bfs3d.rea
> -3 spatial dimension (if negative dump .re2 file)
> 1 number of fields
> #
> #=======================================================================
> #
> box_1
> -18 -2 -4 nelx,nely,nelz (negative --> auto spacing)
> 0.0 20.0 1.085 x_0 x_1 ratio
> -1.0 -0.5 2.6
> 0.0 4.0 1
> W ,O ,W ,E ,P ,P bc's ! west,east,south,north,bottom,top
> box_2
> -18 -2 -4 nelx,nely,nelz (negative --> auto spacing)
> 0.0 20.0 1.085 x_0 x_1 ratio
> -0.5 0.0 0.39
> 0.0 4.0 1
> W ,O ,E ,E ,P ,P bc's ! west,east,south,north,bottom,top
> box_3
> -18 -4 -4 nelx,nely,nelz (negative --> auto spacing)
> 0.0 20.0 1.085 x_0 x_1 ratio
> 0.0 5.0 2.8
> 0.0 4.0 1
> E ,O ,E ,SYM,P ,P bc's ! west,east,south,north,bottom,top
> box_4
> -5 -4 -4 nelx,nely,nelz (negative --> auto spacing)
> -10.0 0.0 0.58 x_0 x_1 ratio
> 0.0 5.0 2.8
> 0.0 4.0 1
> v ,E ,W ,SYM,P ,P bc's ! west,east,south,north,bottom,top
> *************************************************************************
>
>
>
> *************************************************************************
> c-----------------------------------------------------------------------
> subroutine planar_average_t(ua,u,w1,w2)
> c
> c Compute r-s planar average of quantity u()
> c
> include 'SIZE'
> include 'GEOM'
> include 'PARALLEL'
> include 'WZ'
> include 'ZPER'
>
> real ua(nz1,nelz),u(nx1,ny1,nx1,nelv),w1(nz1,nelz),w2(nz1,nelz)
> integer e,eg,ex,ey,ez
>
> nz = nz1*nelz
> call rzero(ua,nz)
> call rzero(w1,nz)
>
> blk1 = 144 ! 18*2*4
> blk2 = 288 ! 18*2*4
> blk3 = 576 ! 18*4*4
> blk4 = 656 ! 18*4*4
>
> do e=1,nelt
> eg = lglel(e)
> if(eg .le. blk1) then
> nelx = 18
> nely = 2
> nelz = 4
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
> elseif(eg .gt. blk1 .and. eg .le. blk2) then
> nelx = 18
> nely = 2
> nelz = 4
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
> elseif(eg .gt. blk2 .and. eg .le. blk3) then
> nelx = 18
> nely = 4
> nelz = 4
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
> elseif(eg .gt. blk3 .and. eg .le. blk4) then
> nelx = 5
> nely = 4
> nelz = 4
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
> endif
>
> do k=1,nz1
> do j=1,ny1
> do i=1,nx1
> zz = (1.-zgm1(k,3))/2. ! = 1 for i=1, = 0 for k=nx1
> aa = zz*area(i,1,5,e) + (1-zz)*area(i,1,6,e) ! wgtd jacobian
> w1(k,ez) = w1(k,ez) + aa
> ua(k,ez) = ua(k,ez) + aa*u(i,j,k,e)
> enddo
> enddo
> enddo
> enddo
> call gop(ua,w2,'+ ',nz)
> call gop(ua,w2,'+ ',nz)
>
> do i=1,nz
> ua(i,1) = u(a,i)/w1(i,1) ! Normalize
> enddo
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine planar_fill_t(u,ua)
> c
> c Fill array u with planar values from ua().
> c For tensor-product array of spectral elements
> c
> include 'SIZE'
> include 'GEOM'
> include 'PARALLEL'
> include 'WZ'
> include 'ZPER'
>
> real u(nx1,ny1,nz1,nelv),ua(lz1,lelz)
> integer e,eg,ex,ey,ez
>
> blk1 = 144 ! 18*2*4
> blk2 = 288 ! 18*2*4
> blk3 = 576 ! 18*4*4
> blk4 = 656 ! 18*4*4
>
> do e=1,nelt
> eg = lglel(e)
> if(eg .le. blk1) then
> nelx = 18
> nely = 2
> nelz = 4
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
> elseif(eg .gt. blk1 .and. eg .le. blk2) then
> nelx = 18
> nely = 2
> nelz = 4
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
> elseif(eg .gt. blk2 .and. eg .le. blk3) then
> nelx = 18
> nely = 4
> nelz = 4
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
> elseif(eg .gt. blk3 .and. eg .le. blk4) then
> nelx = 5
> nely = 4
> nelz = 4
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
> endif
>
> do k=1,nz1
> do j=1,ny1
> do i=1,nx1
> u(i,j,k,e) = ua(k,ez)
> enddo
> enddo
> enddo
> enddo
>
> return
> end
> c-----------------------------------------------------------------------
> *************************************************************************
>
>
More information about the Nek5000-users
mailing list