[Nek5000-users] middle plane flux

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

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

----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.

-------------- next part --------------
c-----------------------------------------------------------------------
      subroutine uservp (ix,iy,iz,ieg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
C
      udiff =0.
      utrans=0.
      return
      end
c-----------------------------------------------------------------------
      subroutine userf  (ix,iy,iz,ieg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
c      include 'TRIPF'

      ffx = 0.0 
      ffy = 0.0
      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    ! called once per step
      include 'SIZE'
      include 'TOTAL'
       if (istep.gt.10) call flux_2d(vz)

      return
      end
c-----------------------------------------------------------------------
      subroutine userbc (ix,iy,iz,iside,ieg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

 

      ux=0.0
      uy=0.0
      uz=0.0

      return
      end
c-----------------------------------------------------------------------
      subroutine useric (ix,iy,iz,ieg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      real pi,zs
      real a,s

      ux = 0.
      uy = 0.
      uz = 3.0*(2.0-x**2-y**2)/4.0

      s = 1.0
     
      
      call random_number(a)
      a=a-0.5
      ux=ux+s*a
      
      call random_number(a)
      a=a-0.5
      uy=uy+s*a
      
      call random_number(a)
      a=a-0.5
      uz=uz+s*a
      
      
      return
      end

C=======================================================================


      SUBROUTINE usrdat
      include 'SIZE'
      include 'TOTAL'
      include 'ZPER'    ! For nelx,nely,nelz - needed for z_average


      nelx = 22
      nely = 22
      nelz = 95

      RETURN
      END

C=======================================================================

      subroutine usrdat2
      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat3
      return
      end
c----------------------------------------------------------------------
c-----------------------------------------------------------------------
      subroutine flux_2d(uz)

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

      real uz(1)
      common /scravg/ w(1),w1(1),w2(1)

      real Re_bulk, Re_mid
      Re_bulk = 1.0 / param(2)


      call zy_plane_average(w,uz,w1,w2)

      Re_mid = Re_bulk * w(1)

      if (nid.eq.0) then	 	
         open(unit=37,access='append',file='history')
c	 write(37,'(A)') '  Time   <flux>   '
	 write(37,3) time,Re_mid     
         if (nid.eq.0.and.istep.eq.NTSTEPS) then
            close(37)
         endif
         
 3       format(1p15e17.9)
      endif
c      write(*,*) 'mid flux:::::::::', w

      return
      end
c-----------------------------------------------------------------------
      subroutine zy_plane_average(ua,u,w1,w2)
c
c     Compute the zy average of quantity u() in the mid plane of the domain in the x direction.
c
      include 'SIZE'
      include 'GEOM'
      include 'PARALLEL'
      include 'WZ'
      include 'ZPER'

      real ua(1),u (nx1,ny1,nz1,nelv)
     $    ,w1(1),w2(1)

      integer e,eg,ex,ey,ez
      real dy2,dz2
      integer frame, plane

C     Find the mid-plane element-number for even/odd number of elements in the x direction
      if (mod(nelx,2).eq.0) then
         frame = (real(nelx) / 2.0) + 1
         plane = 1
      else if (mod(nelx,2).eq.1) then
         frame = floor(real(nelx) / 2.0)
         plane = floor(real(nx1) / 2.0) + 1
      end if
 
      nelxy = nelx*nely
      if (nelxy.gt.lelx*lely) call exitti
     $ ('ABORT IN yz_plane_average.Increase lelx*lely in SIZE:$',nelxy)

      call rzero(ua,1)
      call rzero(w1,1)

      do e=1,nelt

         eg = lglel(e)

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

         i = plane
         if (ex.eq.frame) then
         do j=1,ny1
         do i=plane,plane
            dy2 = 1.0*( ym1(i,ny1,k,e) - ym1(i,1,k,e) ) ! Computing the length of the elements in the y direction
            dz2 = 1.0                                   ! Assuming uniform in "z" direction
            do k=1,nz1

               ua(i)=ua(i)+dy2*wym1(j)*dz2*wzm1(k)*u(i,j,k,e)
               w1(i)=w1(i)+dy2*wym1(j)*dz2*wzm1(k) ! the weigths which must be equal to the yz area


            enddo
         enddo
         enddo
         endif
      enddo

      call gop(ua,w2,'+  ',mxy)
      call gop(w1,w2,'+  ',mxy)

         ua(:) = ua(:) / w1(:)   ! Normalize 

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

c automatically added by makenek
      subroutine usrsetvert(glo_num,nel,nx,ny,nz) ! to modify glo_num
      integer*8 glo_num(1)
      return
      end