[Nek5000-users] Area integral calculation in Nek

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

Thanks Joerg,

Looks good.  One thing to note that nz1 is variable and not a parameter,
so there will be an issue with the fbar() declaration.
I would suggest using lz1 and lelz for declaration of fbar(),
etc., as these parameters are already established for that role.
I think the issue you ran into before was that I failed to include
ZPER in the include statement, which is how nelz becomes known
to the routine.  (It remains, however, up to the user to set
nelz and lelz to the proper values.)

Also, I often find it handy to define a "zbar()" which computes
the mean of z and then prints this out along with fbar and wght.
That makes it easy to plot (zbar,fbar), etc.

Cheers,

Paul


On Thu, 18 Oct 2012, nek5000-users at lists.mcs.anl.gov wrote:

> 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
>
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>