[Nek5000-users] middle plane flux

[nek5000-users at lists.mcs.anl.gov]

Dear NEKs

I think I could fix the problem regarding mid-plane flux computation
somehow. I have attached the routine in case any users run into the same
issue in the future!

Thanks
Azad


P.S. Just wanted to mention I add the plane number manually (if
(ex.eq.4) then) since I had 6 element in the x-direction. One could
change it with respect to the domain.


> Date: Thu, 01 Mar 2012 22:09:57 +0100
> From: nek5000-users at lists.mcs.anl.gov
> Subject: [Nek5000-users] middle plane flux
> To: nek5000-users at lists.mcs.anl.gov
> Message-ID: <20120301220957.994227o61in5hbj9 at www.mech.kth.se>
> Content-Type: text/plain; charset="iso-8859-1"; Format="flowed";
> 	DelSp="Yes"
> 
> Dear NEKs
> 
> I am using NEK5000 to simulate turbulent flow in duct (channel with  
> side walls). I would like to calculate flux at the middle slab of the  
> domain (a zy plane) in each time step in order to compare it with my  
> bulk (Z is my streamwise direction). Attached file is the .usr file  
> containing the routine I wrote to compute the mid-flux. Seems works  
> fine for the primitive tests; however, when I increase number of  
> processors to 256 or upper it gives me NaNs. I was wondering if you  
> could help me with that.
> 
> With many thanks in advance.
> Azad

> 
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
> 
> 
> End of Nek5000-users Digest, Vol 37, Issue 1
> ********************************************

-------------- next part --------------
c-----------------------------------------------------------------------
      subroutine userchk    
      include 'SIZE'
      include 'TOTAL'

      call midplane_flux(vx,vy,vz)

      return
      end
c-----------------------------------------------------------------------
      subroutine midplane_flux(ux,uy,uz)

      include 'SIZE'
      include 'TOTAL'
      include 'ZPER'

      real ux(1),uy(1),uz(1)

      parameter (xsli = ly1*lz1*lely*lelz)
      parameter (yzaver = 1)
      common /x_sli/ u(xsli),v(xsli),w(xsli),w1(xsli),w2(xsli)
      common /yz_aveg/uuu(yzaver),vvv(yzaver),www(yzaver),
     $     www1(yzaver),www2(yzaver)


      call my_x_slice(w,uz,w1,w2)
      call my_yz_average(www,w,www1,www2)



      return
      end
c-----------------------------------------------------------------------
      subroutine my_yz_average(ua,u,www1,www2)

      include 'SIZE'
      include 'GEOM'
      include 'PARALLEL'
      include 'WZ'
      include 'ZPER'

      real ua(1),u (ny1,nz1,nely*nelz) 
     $    ,www1(1),www2(1)
      integer e,eg,ex,ey,ez
      real dz2,dy2


      mxy = 1
      call rzero(ua,mxy)
      call rzero(www1,mxy)

      do e=1,nely*nelz

         eg = lglel(e)
         call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)


         do j=1,ny1
            dy2 = 1.0*( ym1(1,ny1,1,e) - ym1(1,1,1,e) )
            dz2 = 1.0  !  Assuming uniform in "z" direction
            do k=1,nz1
               ua(1) = ua(1)+dy2*wym1(j)*dz2*wzm1(k)*u(j,k,e)
               www1(1) = www1(1)+dy2*wym1(j)*dz2*wzm1(k) 
            enddo
         enddo
      enddo

      call gop(ua,www2,'+  ',mxy)
      call gop(www1,www2,'+  ',mxy)

      do i=1,mxy
         ua(i) = ua(i) / www1(i)   ! Normalize
      enddo

      return
      end

c-----------------------------------------------------------------------
      subroutine my_x_slice (ua,u,w1,w2)

      include 'SIZE'
      include 'GEOM'
      include 'PARALLEL'
      include 'WZ'
      include 'ZPER'
c
      real ua(ny1,nz1,nely,nelz),u (nx1,ny1,nz1,nelv)
     $    ,w1(ny1,nz1,nely,nelz),w2(ny1,nz1,nely,nelz)
      integer e,eg,ex,ey,ez
      real dx2
c
      myz = nely*nelz*ny1*nz1
      call rzero(ua,myz)
c
      do e=1,nelt
c
         eg = lglel(e)
         call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)

         i = 1
         if (ex.eq.4) then
            do k=1,nz1
            do j=1,ny1
               ua(j,k,ey,ez) = u(i,j,k,e)
            enddo
            enddo
         endif
      enddo

      call gop(ua,w2,'+  ',myz)

      return
      end
c-----------------------------------------------------------------------