[Nek5000-users] Two-way coupled LPT

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

Hi Neks,


I'm currently working on improving my two-way coupling set up for my Lagrangian Particle Tracker which I'm developing as part of my PhD.

To do this I require to know the coordinates of the corners of each cell (the space around a mesh gridpoint), in order to determine whether a particle is contained in that cell.


I am having trouble doing this for the entire mesh, and I think it's something to do with the parallel processing. My approach is as follows:


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


if(TWOWAYON.EQ.1) then
          do e = 1, nelv
            ieg = lglel(e)
            write(*,*) 'IEG:',e
            write(*,*) 'nelv:',nelv
            write(*,*) 'nelt:',nelt
            do i = 1, nx1*nelx
              do j = 1, ny1*nely
                do k = 1, nz1*nelz
c                  write(*,*) 'X:',i,'Y:',j,'Z:',k

                   if(i.GE.2) then
                     xmc(i,j,k,e) = xm1(i,j,k,e) -
     $               (xm1(i,j,k,e)-xm1(i-1,j,k,e))*0.5D0
                   else
                     xmc(i,j,k,e)  = xm1(i,j,k,e)
                   endif

                   if(i.LE.7) then
                     xxc(i,j,k,e) = xm1(i,j,k,e) +
     $               (xm1(i+1,j,k,e)-xm1(i,j,k,e))*0.5D0
                   else
                     xxc(i,j,k,e) = xm1(i,j,k,e)
                   endif

                   if(j.GE.2) then
                     ymc(i,j,k,e) = ym1(i,j,k,e) -
     $               (ym1(i,j,k,e)-ym1(i,j-1,k,e))*0.5D0
                   else
                     ymc(i,j,k,e) = ym1(i,j,k,e)
                   endif

                   if(j.LE.7) then
                     yxc(i,j,k,e) = ym1(i,j,k,e) +
     $               (ym1(i,j+1,k,e)-ym1(i,j,k,e))*0.5D0
                   else
                     yxc(i,j,k,e) = ym1(i,j,k,e)
                   endif

                   if(k.GE.2) then
                     zmc(i,j,k,e) = zm1(i,j,k,e) -
     $               (zm1(i,j,k,e)-zm1(i,j,k-1,e))*0.5D0
                   else
                     zmc(i,j,k,e) = zm1(i,j,k,e)
                   endif

                   if(k.LE.7) then
                     zxc(i,j,k,e) = zm1(i,j,k,e) +
     $               (zm1(i,j,k+1,e)-zm1(i,j,k,e))*0.5D0
                   else
                     zxc(i,j,k,e) = zm1(i,j,k,e)
                   endif

                   vol(i,j,k,e) =
     $             (xxc(i,j,k,e) - xmc(i,j,k,e))
     $            *(yxc(i,j,k,e) - ymc(i,j,k,e))
     $            *(zxc(i,j,k,e) - zmc(i,j,k,e))

c                  write(*,*) 'Xmin: ', xmc(i,j,k,e)
c                  write(*,*) 'Xmax: ', xxc(i,j,k,e)
c                  write(*,*) 'Ymin: ', ymc(i,j,k,e)
c                  write(*,*) 'Ymax: ', yxc(i,j,k,e)
c                  write(*,*) 'Zmin: ', zmc(i,j,k,e)
c                  write(*,*) 'Zmax: ', zxc(i,j,k,e)
c                  write(*,*) 'VOL : ', vol(i,j,k,e)
                enddo
              enddo
            enddo
          enddo
        endif

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


Now if I sum up the volume using: tvol = glsum(vol,n), I get a weird volume that only seems to represent one element (ish).


If you can see where I'm going wrong I'd greatly appreciate it.


Thanks,


Lee Mortimer

University of Leeds
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20160129/e4776c6e/attachment.html>