[Nek5000-users] Probability Density Function

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

Hi Joerg,

For a weight, I would just use the diagonal mass matrix, bm1(i,j,k,e) ,
which corresponds precisely to the volume associated with a gridpoint.

I don't know why your weights went negative --- assuming that your
mesh was generated with n2to3, it would seem that they should indeed
be positive.


Paul


On Mon, 17 Dec 2012, nek5000-users at lists.mcs.anl.gov wrote:

>
>
> Dear Neks,
>
> I want to calculate the probability density function of
> a scalar field ff in the attached subroutine.
> Since the element mesh is
> nonuniform I have to assign a weight "wght" with each grid point of
> ff.
> The weight is given by the volume of the associated grid cell.
> Coordinates of x,y,z of all vertices
> (primary and secondary nodes) are
> stored in xm1,ym1,zm1 as far as I know.
> This weight turns out to have
> in parts negative values which is unphysical. What went wrong?
>
>
> subroutine pdf_calc(ff,step,i_offset,i_name)
> include 'SIZE'
> include
> 'TOTAL'
>
> parameter(npdf=501)
>
> real ff(lx1,ly1,lz1,lelt)
> real
> pdf(npdf)
> real work(npdf)
> real val, vol, offset, wght
>
> integer
> e,eg,ex,ey,ez,f
>
> !-----Set arrays to zero
> call rzero(pdf,npdf)
> call
> rzero(work,npdf)
>
> !-----Offset
> if(i_offset==0) offset=0.0
>
> if(i_offset==1) offset=int(npdf/2)*step
>
> !-----Pick face 5
> f = 5
>
>
> vol=atan(1.0)
>
> do e=1,nelv
> do k=1,nz1
> do j=1,ny1
> do i=1,nx1
> wght=
> (zm1(i,j,k+1,e)-zm1(i,j,k,e)) * area(i,j,f,e)  mpi_allreduce)
> call
> gop(pdf,work,'+
> ',npdf)
>
> !------------------------------------------------------------
>
> if(nid.eq.0)then
>
>
> if(i_name.eq.1)OPEN(10,file="pdf_uzte.dat",position="append")
>
> if(i_name.eq.2)OPEN(10,file="pdf_epst.dat",position="append")
>
> if(i_name.eq.3)OPEN(10,file="pdf_epsv.dat",position="append")
>
> if(i_name.eq.4)OPEN(10,file="pdf_temp.dat",position="append")
>
> if(i_name.eq.5)OPEN(10,file="pdf_dtdz.dat",position="append")
>
> do
> ipdf=1,npdf
> write(10,*) (ipdf-1)*step-offset, pdf(ipdf)
> enddo
>
> CLOSE(10)
>
> endif
>
> return
> end
>
> Thanks in advance and best wishes,
> Joerg.
>
> ----------------------------------------------------
> Joerg
> Schumacher
>
> Heisenberg Professor for Theoretical Fluid
> Mechanics
> Institute of Thermodynamics and Fluid Mechanics
> Department of
> Mechanical Engineering
> Ilmenau University of Technology
> P.O. Box 100565
>
> D-98684 Ilmenau
> Germany
>
> E-mail:
> joerg.schumacher at tu-ilmenau.de
> http://www.tu-ilmenau.de/tsm
> Phone:
> +49-3677-69-2428
> Fax: +49-3677-69-2411
>