[Nek5000-users] Area integral calculation in Nek

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

Hi Paul,

I attach a version with two corrections. This one should work now correctly.
Thanks once more.

Joerg.

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
      subroutine vertical_mean_profile(ff)
      include 'SIZE'
      include 'TOTAL'

      integer nelz
      parameter(nelz=32)

      real ff(nx1,ny1,nz1,nelv)

      real fbar(1:nz1,1:nelz)
      real wght(1:nz1,1:nelz)
      real work(1:nz1,1:nelz)

      integer e,eg,ex,ey,ez,f

      nelxy = nelgv/nelz 

!-----Set both arrays to zero
      call rzero(fbar,nz1*nelz)
      call rzero(wght,nz1*nelz)


!-----Pick face 5 to evaluate surface Jacobian
      f = 5 

      do e=1,nelv

         eg = lglel(e)
         call get_exyz(ex,ey,ez,eg,nelxy,1,nelz)

         do k=1,nz1
           do i=1,nx1*ny1
             fbar(k,ez) = fbar(k,ez)+area(i,1,f,e)*ff(i,1,k,e)
             wght(k,ez) = wght(k,ez)+area(i,1,f,e)
           enddo
         enddo

      enddo

!-----Gather over all processes (-> mpi_allreduce)
      call gop(fbar,work,'+  ',nz1*nelz)
      call gop(wght,work,'+  ',nz1*nelz)

!-----Area average
      do i=1,nz1*nelz
        fbar(i,1)=fbar(i,1)/wght(i,1) 
      enddo

!-----Output of the vertical profile of planar averages 
!
!     Example with 32 (=nelz) vertical elements and lx1=ly1=lz1=4:
!      -> 2 internal vertical GLL points per element
!      -> 97 vertical planes, here output of 96
!
!------------------------------------------------------------ 
      if(nid.eq.0)then
        do k=1,nelz
        do i=1,nz1-1
          write(10,*)fbar(i,k),wght(i,k)
        enddo   
        enddo
      endif

      return
      end




Am 15.10.2012 um 22:15 schrieb nek5000-users at lists.mcs.anl.gov:

> 
> Hi Joerg,
> 
> If your geometry came from n2to3, then your data
> is organized into planar sections and you can compute
> 
>           /
>  F(z) =   |  f(x,y,z) dx dy
>           /A
> 
> as follows:
> 
>      common /scruz/ fbar(lz1,lelz),wght(lz1,lelz),work(lz1,lelz)
>      integer e,eg,ex,ey,ez,f
> 
>      nelz  = [how many levels you have in z]
>      nelxy = nelgv/nelz ! Number of elements in a plane
> 
>      call rzero(fbar,nz1*nelz)
>      call rzero(wght,nz1*nelz)
> 
>      f = 5  ! Pick face 5 to evaluate surface Jacobian
> 
>      do e=1,nelv
>         eg = lglel(e)
>         call get_exyz(ex,ey,ez,eg,nelxyz,1,nelz)
>         do k=1,nz1
>         do i=1,nx1*ny1
>            fbar(k,ez) = fbar(k,ez)+area(i,f,e)*f(i,1,k,e)
>            wght(k,ez) = wght(k,ez)+area(i,f,e)
>         enddo
>         enddo
>      enddo
>      call gop(fbar,work,'+  ',nz1*nelv)
>      call gop(wght,work,'+  ',nz1*nelv)
>      do i=1,nz1*nelz
>         fbar(i,1)=fbar(i,1)/wght(i,1) ! If you want the average
>      enddo
> 
> That should pretty much work as is... You need to specify
> nelz and lelz >= nelz.
> 
> Paul
> 
> 
> On Mon, 15 Oct 2012, nek5000-users at lists.mcs.anl.gov wrote:
> 
>> Hi,
>> 
>> I work in cylindrical geometry and want to evaluate an integral of one of
>> scalar fields over the circular cross section at fixed position z_0. I
>> suppose that the best way would be to make use of the GLL quadrature formulas
>> on the elements after picking the elements at desired height z_0. I have two
>> short questions:
>> 
>> 1/  Is it correct that xm1, ym1, zm1 are the coordinate arrays, i.e. fixing
>> for example zm1(i,1,1,1)=0.5 would then give all elements that contain nodes
>> in exactly this plane?
>> 
>> 2/ Are the GLL weights at the primary and secondary nodes stored in a
>> specific array or do I have to evaluate them for each element, e.g. by
>> combination of 1d routines which are collected in speclib.f?
>> 
>> Thanks in advance, Joerg.
>> 
>> _______________________________________________
>> Nek5000-users mailing list
>> Nek5000-users at lists.mcs.anl.gov
>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>> 
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users