[Nek5000-users] Planar average in z-direction problem (image updated)

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

Dear Hyunduk,


I see several issues with the code you have.


The first point is that your code will not work properly

in parallel, so you need to be careful about parallel treatment

of IO if you are running on more than one MPI rank.


Secondly, it's generally not a good idea to test for equality

with real numbers -- are you certain that the y values of

interest are 0.0000000 ? or are they potentially 1.e-7 (say).

Some of the mesh generation tools will have noise and your

zero values might not be exactly zero.


zgm1 is the position of the nodal points on the interval [-1,1]

in the reference element


area() is the surface area at the top and bottom of each

element.   (Here, we assume that "top" in z means also

top in r-s-t space for the canonical Omega-hat := [-1,1]^3.)


Thus, aa is a linear combination of the area at the top of

the given element and bottom of the given element.   This

way the surface integral can handle tapered elements

(but nothing more complex than that).


hth,


Paul




________________________________
From: Nek5000-users <nek5000-users-bounces at lists.mcs.anl.gov> on behalf of nek5000-users at lists.mcs.anl.gov <nek5000-users at lists.mcs.anl.gov>
Sent: Tuesday, July 24, 2018 1:34 AM
To: nek5000-users at lists.mcs.anl.gov
Subject: [Nek5000-users] Planar average in z-direction problem (image updated)


Hi, Neks



I wonder why planar_average_z in navier5.f shows wrong result.

here I attached 2 images.



1, https://drive.google.com/file/d/1Kng9HZfQW-TGTFwmFcXzbrcvhRYrQSdc/view?usp=sharing

2. https://drive.google.com/file/d/1FN3i5mE6VSoPeehG5evjC33WrD3J4lB5/view?usp=sharing





First image shows the z-directional grids points from zm1

(I wrote a code to extract z coordinate at y = 0 as follows:

======================================================

insnum = nid

insuni = insnum + 1000

WRITE(insnam,98) 'inst_field_', insnum

OPEN(unit=insuni,file=insnam

      insnum = nid

      insuni = insnum + 1000

      WRITE(insnam,98) 'inst_field_', insnum

      Open(unit=insuni,file=insnam)



      do i = 1,n

      if (ym1(i,1,1,1).eq.0.0)then

      WRITE(insuni,5) xm1(i,1,1,1)

     %                  ,  zm1(i,1,1,1)

     %                  ,  vx(i,1,1,1)

     %                  ,  vy(i,1,1,1)

     %                  ,  vz(i,1,1,1)

     %                  ,  phig(i,1,1,1)

      endif



      enddo

      CLOSE(insuni)

======================================================



Then, I can collect lots of field data from each processor and sort data with ascending order using matlab.

So, I got the result as 1st image.



But when I use plannar_average_z

I got the result as shown in 2nd image.

1. I really want to modify the problems, so could you give me your idea on this problem?

2. Additionally, could you explain 'zgm1' and 'area'?

I noticed that, planar_average_r or s ... in turbChannel case has also 'area' but, those case are used to calculate on different direction.
But they used same input parameters.



What I understand about area is, surface area of element.

But I thought that different direction needs different indexing



=======================================================

      subroutine planar_average_z(ua,u,w1,w2)

c

c     Compute r-s planar average of quantity u()

c

      include 'SIZE'

      include 'GEOM'

      include 'PARALLEL'

      include 'WZ'

      include 'ZPER'

c

      real ua(lz1,nelz),u(lx1*ly1,lz1,nelv),w1(lz1,nelz),w2(lz1,nelz)

      integer e,eg,ez

c

      melxy = nelx*nely

c

      nz = lz1*nelz

      call rzero(ua,nz)

      call rzero(w1,nz)

c

      do e=1,nelt

c

          eg = lglel(e)

c         ez = 1 + (eg-1)/melxy

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

c

         do k=1,lz1

         do i=1,lx1*ly1

            zz = (1.-zgm1(k,3))/2.  ! = 1 for k=1, = 0 for k=lz1

            aa = zz*area(i,1,5,e) + (1-zz)*area(i,1,6,e)  ! wgtd jacobian

            w1(k,ez) = w1(k,ez) + aa

            ua(k,ez) = ua(k,ez) + aa*u(i,k,e)

         enddo

         enddo

      enddo

c

      call gop(ua,w2,'+  ',nz)

      call gop(w1,w2,'+  ',nz)

c

      do i=1,nz

         ua(i,1) = ua(i,1) / w1(i,1)   ! Normalize

      enddo

c

      return

      end

======================================================



Regards,



Hyunduk.

------------------------------------------------------------------------------------------------------------------

Hyunduk Seo (Mr.)

Integrated Ph.D. course, School of Mechanical Engineering, Pusan National University

Experimental Thermo-Fluids Mechanics and Energy Systems (ExTENsys) Laboratory

http://fluid.pusan.ac.kr<http://fluid.pusan.ac.kr/>
T.: +82-51-510-1536 M.: +82-10-2964-0994



[http://webmail.pusan.ac.kr/api/notify.php?home=MjAxODA3MjQxNTM0MTR8bmVrNTAwMC11c2Vyc0BsaXN0cy5tY3MuYW5sLmdvdnxoZHNlbzAxMjNAcHVzYW4uYWMua3I=]
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180724/31a6d3b7/attachment.html>