<div name='mp-default' style='font-family:돋움,Dotum;font-size:10pt;'><p>Dear Paul,</p><p> </p><p>Thank you for you recommendation. </p><p> </p><p>As you said, I also understand that my code is not MPI friendly.</p><p> </p><p>But, my topic was malfunction of the planar_average_z in my case and I just wanted to know whether there are any problems in my grid points.</p><p> </p><p>In fact, I noticed that my computational domains with 4 boxes made from .box file have some problems.</p><p>(So, I wrote a code to generate mesh grid for 1 box with clustering. And I manually assigned grid points then, problem was solved.)</p><p> </p><p>I think connectivity between each box was bad at the time.</p><p> </p><p>Regards,</p><p> </p><p>Hyunduk.</p><p> </p><p>---------------------------------------------------------------------------</p><p>---------------------------------------------------------------------------</p><p>Title: Re: Planar average in z-direction problem (image updated)</p><p>Dear Hyunduk,</p><p> </p><p> </p><p>I see several issues with the code you have.</p><p> </p><p> </p><p>The first point is that your code will not work properly</p><p> </p><p>in parallel, so you need to be careful about parallel treatment</p><p> </p><p>of IO if you are running on more than one MPI rank.</p><p> </p><p> </p><p>Secondly, it's generally not a good idea to test for equality</p><p> </p><p>with real numbers -- are you certain that the y values of</p><p> </p><p>interest are 0.0000000 ? or are they potentially 1.e-7 (say).</p><p> </p><p>Some of the mesh generation tools will have noise and your</p><p> </p><p>zero values might not be exactly zero.</p><p> </p><p> </p><p>zgm1 is the position of the nodal points on the interval [-1,1]</p><p> </p><p>in the reference element</p><p> </p><p> </p><p>area() is the surface area at the top and bottom of each</p><p> </p><p>element.   (Here, we assume that "top" in z means also</p><p> </p><p>top in r-s-t space for the canonical Omega-hat := [-1,1]^3.)</p><p> </p><p> </p><p>Thus, aa is a linear combination of the area at the top of</p><p> </p><p>the given element and bottom of the given element.   This</p><p> </p><p>way the surface integral can handle tapered elements</p><p> </p><p>(but nothing more complex than that).</p><p> </p><p> </p><p>hth,</p><p> </p><p> </p><p>Paul</p><p> </p><p>---------------------------------------------------------------------------</p><p>---------------------------------------------------------------------------</p><p> </p><p>Title: Planar average in z-direction problem (image updated)</p><p> </p><p>Hi, Neks</p><p> </p><p> </p><p> </p><p>I wonder why planar_average_z in navier5.f shows wrong result.</p><p> </p><p>here I attached 2 images.</p><p> </p><p> </p><p> </p><p>1, https://drive.google.com/file/d/1Kng9HZfQW-TGTFwmFcXzbrcvhRYrQSdc/view?usp=sharing</p><p> </p><p>2. https://drive.google.com/file/d/1FN3i5mE6VSoPeehG5evjC33WrD3J4lB5/view?usp=sharing</p><p> </p><p> </p><p> </p><p> </p><p> </p><p>First image shows the z-directional grids points from zm1</p><p> </p><p>(I wrote a code to extract z coordinate at y = 0 as follows:</p><p> </p><p>======================================================</p><p> </p><p>insnum = nid</p><p> </p><p>insuni = insnum + 1000</p><p> </p><p>WRITE(insnam,98) 'inst_field_', insnum</p><p> </p><p>OPEN(unit=insuni,file=insnam</p><p> </p><p>      insnum = nid</p><p> </p><p>      insuni = insnum + 1000</p><p> </p><p>      WRITE(insnam,98) 'inst_field_', insnum</p><p> </p><p>      Open(unit=insuni,file=insnam)</p><p> </p><p> </p><p> </p><p>      do i = 1,n</p><p> </p><p>      if (ym1(i,1,1,1).eq.0.0)then</p><p> </p><p>      WRITE(insuni,5) xm1(i,1,1,1)</p><p> </p><p>     %                  ,  zm1(i,1,1,1)</p><p> </p><p>     %                  ,  vx(i,1,1,1)</p><p> </p><p>     %                  ,  vy(i,1,1,1)</p><p> </p><p>     %                  ,  vz(i,1,1,1)</p><p> </p><p>     %                  ,  phig(i,1,1,1)</p><p> </p><p>      endif</p><p> </p><p> </p><p> </p><p>      enddo</p><p> </p><p>      CLOSE(insuni)</p><p> </p><p>======================================================</p><p> </p><p> </p><p> </p><p>Then, I can collect lots of field data from each processor and sort data with ascending order using matlab.</p><p> </p><p>So, I got the result as 1st image.</p><p> </p><p> </p><p> </p><p>But when I use plannar_average_z</p><p> </p><p>I got the result as shown in 2nd image.</p><p> </p><p>1. I really want to modify the problems, so could you give me your idea on this problem?</p><p> </p><p>2. Additionally, could you explain 'zgm1' and 'area'?</p><p> </p><p>I noticed that, planar_average_r or s ... in turbChannel case has also 'area' but, those case are used to calculate on different direction.</p><p>But they used same input parameters.</p><p> </p><p> </p><p> </p><p>What I understand about area is, surface area of element.</p><p> </p><p>But I thought that different direction needs different indexing</p><p> </p><p> </p><p> </p><p>=======================================================</p><p> </p><p>      subroutine planar_average_z(ua,u,w1,w2)</p><p> </p><p>c</p><p> </p><p>c     Compute r-s planar average of quantity u()</p><p> </p><p>c</p><p> </p><p>      include 'SIZE'</p><p> </p><p>      include 'GEOM'</p><p> </p><p>      include 'PARALLEL'</p><p> </p><p>      include 'WZ'</p><p> </p><p>      include 'ZPER'</p><p> </p><p>c</p><p> </p><p>      real ua(lz1,nelz),u(lx1*ly1,lz1,nelv),w1(lz1,nelz),w2(lz1,nelz)</p><p> </p><p>      integer e,eg,ez</p><p> </p><p>c</p><p> </p><p>      melxy = nelx*nely</p><p> </p><p>c</p><p> </p><p>      nz = lz1*nelz</p><p> </p><p>      call rzero(ua,nz)</p><p> </p><p>      call rzero(w1,nz)</p><p> </p><p>c</p><p> </p><p>      do e=1,nelt</p><p> </p><p>c</p><p> </p><p>          eg = lglel(e)</p><p> </p><p>c         ez = 1 + (eg-1)/melxy</p><p> </p><p>         call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)</p><p> </p><p>c</p><p> </p><p>         do k=1,lz1</p><p> </p><p>         do i=1,lx1*ly1</p><p> </p><p>            zz = (1.-zgm1(k,3))/2.  ! = 1 for k=1, = 0 for k=lz1</p><p> </p><p>            aa = zz*area(i,1,5,e) + (1-zz)*area(i,1,6,e)  ! wgtd jacobian</p><p> </p><p>            w1(k,ez) = w1(k,ez) + aa</p><p> </p><p>            ua(k,ez) = ua(k,ez) + aa*u(i,k,e)</p><p> </p><p>         enddo</p><p> </p><p>         enddo</p><p> </p><p>      enddo</p><p> </p><p>c</p><p> </p><p>      call gop(ua,w2,'+  ',nz)</p><p> </p><p>      call gop(w1,w2,'+  ',nz)</p><p> </p><p>c</p><p> </p><p>      do i=1,nz</p><p> </p><p>         ua(i,1) = ua(i,1) / w1(i,1)   ! Normalize</p><p> </p><p>      enddo</p><p> </p><p>c</p><p> </p><p>      return</p><p> </p><p>      end</p><p> </p><p>======================================================</p><p> </p><p> </p><p> </p><p>Regards,</p><p> </p><p> </p><p> </p><p>Hyunduk.</p></div><img src="http://webmail.pusan.ac.kr/api/notify.php?home=MjAxODA3MjYyMjA3Mzl8bmVrNTAwMC11c2Vyc0BsaXN0cy5tY3MuYW5sLmdvdnxoZHNlbzAxMjNAcHVzYW4uYWMua3I=">