Lagrangian particles -- Nearest Neighbor (GPEXCHTYPE_NN) interface

[JavierSierraAusin]

The module for exchanging the data associated to particles GPEXCHTYPE_NN is not operational.
Functions like GPart_InitUserSeed use this%gpcomm_%VDBSynch even when GPEXCHTYPE_NN is active.
In a (simple) first attempt to use the nearest neighbours interface I allow the allocation of memory of

         ALLOCATE (this%gptmp0_(3, this%maxparts_))
         ALLOCATE (this%vdb_(3, this%maxparts_))

even when GPEXCHTYPE_NN is active. However, in this case GPartComm_PartExchangeV does not work as it should.
When there isn't any exchange of particles, it adds more particles. That I think it could be solved by adding

    this%rtbuff_ = GPNULL
    this%rbbuff_ = GPNULL
    this%sbbuff_ = GPNULL
    this%stbuff_ = GPNULL

before the comment ! Post receives:.
Another issue is on line 1195, which packs on the variable this%sbbuff_ but it then sends this%stbuff_. Then that line I think should be
CALL GPartComm_PPackV(this,this%stbuff_,this%nbuff_,id,px,py,pz,nparts,this%itop_,this%ntop_)
However, even after these changes, the interface does not work properly. After the exchange of particles between neighbouring nodes, I observe that two particles may share the same global id, which is forbidden, and I also think that GPartComm_ConcatV does not work properly because it does not remove all the particles that have been sent away.

Could you take a look into this?

[pmininni]

Yes! We noticed that there are some issues with GPEXCHTYPE_NN, and we should fix this in the next few days. This will make Lagrangian particles work with GPEXCHTYPE_NN, and we will later fix all other particles. In the meantime, note that for small runs the GPEXCHTYPE_VDB method should be faster (a 1024^3 run with 1 million particles has ~1000 less particles than grid points).

[pmininni]

We just uploaded an updated version that should do GPEXCHTYPE_NN for all types of particles. Please tell us if you find any problem.

[JavierSierraAusin]

Thanks! I am going to test it.

If you allow me, I have another question. In the spline interpolation, when computing in parallel (MPI or hybrid)
in the subroutine GPSplineInt_PartUpdate3D, there is a line this%klg_(1,j) = max(min(this%klg_(1,j),kmax),kmin).

Such a line basically decenter the spline interpolation. If we assume that we are in a cell that is in the position (xp,yp,zp) within a slab of limits 1,nz and zp > nz-1 instead of being around the cell i with stencil (i-1,i,i+1,i+2) we would have a stencil (i-2,i-1,i,i+1). The reason behind this I suspect is the fact that the supposedly extended field this%esplfld_ is not extended and the position i+2 = nz+1 is outsite of the slab. I have tested this and in fact when I remove the line this%klg_(1,j) = max(min(this%klg_(1,j),kmax),kmin) on GPSplineInt_PartUpdate3D the coefficients this%esplfld_ associated to the cell i+2 is "empty". This implies that the following line does not work as it should
CALL this%gpcomm_%SlabDataExchangeSF(this%esplfld_,field).
Having a decentered stencil wouldn't be a problem if the coefficients are correctly computed accounting for that but I think it is not the case. Could you check it?

Thanks again!

[JavierSierraAusin]

Follow up of my message.

The communicator SlabDataExchangeSF(this%esplfld_,field) exchanges 2 ghost positions.
So in the this%myrank_ .EQ. 0 corresponds to indices in the regular grid 1:nz, with positions 0:nz-1; in the extended grid indices 1:nz+4 and positions -2:nz+1. The communicator then exchanges the indices (on the extended grid), 1,2and nz+3,nz+4.

If my previous statement is correct, then the line this%klg_(1,j) = (zp(j)-this%xbnds_(3,1))*this%dxi_(3) should be replaced by this%klg_(1,j) = (zp(j)-this%xbnds_(3,1))*this%dxi_(3) - 1, here this%xbnds_(3,1) = -3 and this%zrk_ (j) = (zp(j)-this%xbnds_(3,1)-1)*this%dxi_(3) - real(this%klg_(1,j),kind=GP). Or to modify nzghost = 2. This should be done, in gparts_mod.f90, from

      CALL this%intop_%GPSplineInt_ctor(3, this%nd_, this%libnds_, this%lxbnds_, &
                                        this%tibnds_, this%intorder_ , this%maxparts_, this%gpcomm_, &
                                        this%htimers_(GPTIME_DATAEX), this%htimers_(GPTIME_TRANSP))

to

      CALL this%intop_%GPSplineInt_ctor(3, this%nd_, this%libnds_, this%lxbnds_, &
                                        this%tibnds_, this%intorder_ - 1, this%maxparts_, this%gpcomm_, &
                                        this%htimers_(GPTIME_DATAEX), this%htimers_(GPTIME_TRANSP))

This would follow the logic of the GPartcomm_ctor,

      CALL this%gpcomm_%GPartComm_ctor(GPCOMM_INTRFC_SF, this%maxparts_, &
                                       this%nd_, this%intorder_ - 1, this%comm_, this%htimers_(GPTIME_COMM))

In this case, if the position is, for instance zp = 0.345, the index in the regular grid is 1and in the extended grid should be 3, so the stencil should be (on the extended grid) (2,3,4,5). Similarly, if zp = (nz-1)+0.345 with this choice the stencil would be in the extended grid of indices (nz+1,nz+2,nz+3,nz+4).

Summing up, if this reasoning is correct. The Lagrangian stencil on the z-direction was always decentered (i,i+1,i+2,i+3) instead of (i-1,i,i+1,i+2).

[JavierSierraAusin]

Btw, I did some tests on the implementation that you did. The results coincide with my implementation, so I think the NN module is now working! Thanks!