[Nek5000-users] Area integral calculation in Nek
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Thu Oct 18 14:15:21 CDT 2012
Hi Paul,
> 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.
change to lz1 is fine.
I run in trouble whenever I use lelz.
The attached version works.
subroutine vertical_mean(ff)
include 'SIZE'
include 'TOTAL'
! include 'ZPER'
integer nelz
parameter(nelz=32)
real ff(lx1,ly1,lz1,lelt)
! common /scruz/ fbar(lz1,lelz),wght(lz1,lelz),work(lz1,lelz)
real fbar(lz1,nelz)
real wght(lz1,nelz)
real work(lz1,nelz)
integer e,eg,ex,ey,ez,f
! lelz=32
! nelz=32
nelxy = nelgv/nelz
!-----Set both arrays to zero
call rzero(fbar,lz1*nelz)
call rzero(wght,lz1*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,'+ ',lz1*nelz)
call gop(wght,work,'+ ',lz1*nelz)
!-----Area average
do i=1,lz1*nelz
fbar(i,1)=fbar(i,1)/wght(i,1)
enddo
!-----Output of the vertical profile of plane 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
OPEN(10,file="uzt.dat",access="append")
do k=1,nelz
do i=1,lz1-1
write(10,*)fbar(i,k),wght(i,k)
enddo
enddo
CLOSE(10)
endif
return
end
Thanks Joerg.
> 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
>>
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
More information about the Nek5000-users
mailing list