[Nek5000-users] Planar averaging in 2D

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

Hi David,

The attached routine was written for planar averaging in y --

Have a look - it should in principle work in 2D as well.
Note that it's not been tested in 2D or 3D as far as I know.
I wrote it for someone a week or two ago (with the same
caveat).  Am posting to the user list in case someone knows
if this is working or has been tested/fixed.

Paul

On Fri, 8 Oct 2010, Aleksandr Obabko wrote:

> Hi David,
>
> You are correct -- planar_avarge_z has to be modified.  I'm finishing the 
> report today so after it I can make a pass at the routine.
>
> Best,
> Aleks
>
>
> On Fri, 8 Oct 2010, David Goluskin wrote:
>
>> Hi Aleks,
>> 
>> A few months ago, you told me how to get horizontally averaged temperature
>> profiles using the planar_average_z routine.  To get horizontal averages 
>> for
>> a 2D run, I am guessing I need to either modify the planar averaging 
>> routine
>> or change the vertical coordinate form y to z, which I don`t know how to
>> do.  What do you suggest?
>> 
>> Thank you,
>> 
>> David
>> 
>
-------------- next part --------------
c-----------------------------------------------------------------------
      subroutine planar_average_y(ua,u2,u,w1,w2)
c
c     Compute r-s planar average of quantity u() and u^2()
c
      include 'SIZE'
      include 'GEOM'
      include 'PARALLEL'
      include 'WZ'
      include 'ZPER'
c
      real ua(ny1,nely),u2(ny1,nely)   ! Output:  < u >  and < u^2 >
      real u(nx1,ny1,nz1,nelv)         ! Input
      real w1(ny1,nely),w2(ny1,nely)   ! work arrays
      integer e,eg,ex,ey,ez

      ny = ny1*nely
      call rzero(ua,ny)
      call rzero(u2,ny)
      call rzero(w1,ny)

      do e=1,nelt

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

         do k=1,nz1
         do i=1,nx1
         do j=1,ny1
            yy = (1.-zgm1(j,2))/2.  ! = 1 for j=1, = 0 for j=ny1
            aa = yy*area(i,1,1,e) + (1-yy)*area(i,1,3,e)  ! wgtd jacobian, fc 1&3
            w1(j,ey) = w1(j,ey) + aa
            ua(j,ey) = ua(j,ey) + aa*u(i,j,k,e)
            u2(j,ey) = u2(j,ey) + aa*u(i,j,k,e)**2
         enddo
         enddo
         enddo
      enddo

      call gop(ua,w2,'+  ',ny)
      call gop(u2,w2,'+  ',ny)
      call gop(w1,w2,'+  ',ny)

      do i=1,ny
         ua(i,1) = ua(i,1) / w1(i,1)   ! Normalize
         u2(i,1) = u2(i,1) / w1(i,1)   ! Normalize
      enddo

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