Problem in calculation of interaction cells
Dear DualSPHysics developers,
while working with DualSPHysics I found a potential problem in the helper methods that calculates the interaction cells.
This applies only to the GPU version as the CPU version has the identical code but uses no unsigned types.
Please have a look at
The problem are the three lines
I am explaining the principle based on the z coordinate.
nc is of type uint4, which means that nc.z is of type unsigned int.
Because of this the entire calculation is cast to unsigned int, which is problematic in the boundary cases.
If cz is the last cell, the calculation should be
So instead of
Debugging using cuda gdb shows that the following code will be called
Please also note that there are two GPU versions of
Best regards,
Daniel
while working with DualSPHysics I found a potential problem in the helper methods that calculates the interaction cells.
This applies only to the GPU version as the CPU version has the identical code but uses no unsigned types.
Please have a look at
//------------------------------------------------------------------------------
/// Devuelve limites de celdas para interaccion.
/// Returns cell limits for the interaction.
//------------------------------------------------------------------------------
__device__ void KerGetInteractionCells(unsigned rcell
,int hdiv,const uint4 &nc,const int3 &cellzero
,int &cxini,int &cxfin,int &yini,int &yfin,int &zini,int &zfin)
{
//-Obtiene limites de interaccion
//-Obtains interaction limits.
const int cx=PC__Cellx(CTE.cellcode,rcell)-cellzero.x;
const int cy=PC__Celly(CTE.cellcode,rcell)-cellzero.y;
const int cz=PC__Cellz(CTE.cellcode,rcell)-cellzero.z;
//-Codigo para hdiv 1 o 2 pero no cero.
//-Code for hdiv 1 or 2 but not zero.
cxini=cx-min(cx,hdiv);
cxfin=cx+min(nc.x-cx-1,hdiv)+1;
yini=cy-min(cy,hdiv);
yfin=cy+min(nc.y-cy-1,hdiv)+1;
zini=cz-min(cz,hdiv);
zfin=cz+min(nc.z-cz-1,hdiv)+1;
}
The problem are the three lines
zfin=cz+min(nc.z-cz-1,hdiv)+1; for x,y,z.I am explaining the principle based on the z coordinate.
nc is of type uint4, which means that nc.z is of type unsigned int.
Because of this the entire calculation is cast to unsigned int, which is problematic in the boundary cases.
If cz is the last cell, the calculation should be
zfin=cz+min(-1,hdiv)+1;, but due to the unsigned implicit cast it will be zfin=cz+min(UINT_MAX,hdiv)+1;.So instead of
zfin=cz+(-1)+1; // zfin == cz it will be zfin=cz+(hdiv)+1; // zfin > cz .Debugging using cuda gdb shows that the following code will be called
│ __MATH_FUNCTIONS_DECL__ unsigned int min(unsigned int a, int b)
│ {
>│ return umin(a, (unsigned int)b);
│ }
Please also note that there are two GPU versions of
KerGetInteractionCells.Best regards,
Daniel
Comments
zfin=cz+min(int(nc.z)-cz-1,hdiv)+1;Thanks a lot! You are right. This is a bug!
We should use int4 instead of unit4 as we do in the CPU version.
Fortunately the online version will not lead to any error since the number of cells in one axis (nc.x/y/z) is always higher than the cell of the particles.
So that nc.z-cz-1 less than zero is not going to appear (nc.z is always higher than cz).
Indeed there are two functions KerGetInteractionCells().
The first one uses the precomputed cell of each particle so the mentioned problem is not going to appear. The second function KerGetInteractionCells() uses a position and in this case, we can have the problem that the cz computed with that option is higher or equal to nc.z. But in the code we are only using the first function that will not lead to the problem (but we have to solve the problem as you pointed). The second option is there only in case other user/developer wants to compute the interaction with some given points.
Once again Daniel, thanks for sharing this with all
Regards
You are very welcome.
Thank you for the detailed information.
Best regards,
Daniel