Open JavierSierraAusin opened 7 months ago
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).
We just uploaded an updated version that should do GPEXCHTYPE_NN for all types of particles. Please tell us if you find any problem.
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!
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,2
and 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 1
and 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)
.
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!
The module for exchanging the data associated to particles GPEXCHTYPE_NN is not operational. Functions like
GPart_InitUserSeed
usethis%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 ofeven 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 addingbefore the comment
! Post receives:
. Another issue is on line 1195, which packs on the variablethis%sbbuff_
but it then sendsthis%stbuff_
. Then that line I think should beCALL 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 thatGPartComm_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?