When reading the code I found something I could not understand. It seems the code is not in accord with the manual.
the code is in JSPHCPU or JSPHGPU.
//-Shifting correction if(shift && shiftposp1.x!=FLT_MAX){ const float massrhop=massp2/velrhop[p2].w; const bool noshift=(boundp2 && (tshifting==SHIFT_NoBound || (tshifting==SHIFT_NoFixed && CODE_GetType(code[p2])==CODE_TYPE_FIXED))); shiftposp1.x=(noshift? FLT_MAX: shiftposp1.x+massrhop*frx); //-For boundary do not use shifting / Con boundary anula shifting. shiftposp1.y+=massrhop*fry; shiftposp1.z+=massrhop*frz; shiftdetectp1-=massrhop*(drx*frx+dry*fry+drz*frz); }
The shiftdetectp should be the div r as in Eq. 22 in the manual, while the shiftposp1 is the sph gradient operator.
Are the + and - signs opposite from those in the manual？ according to my understanding, the shiftdetectp1 is used to compared with the AFST in Eq.(24) so it should be positive. The shiftposp1 is used to mutiply the D in Eq. (17) or (24), and with for the "-" in in Eq. (17) or (24), the shiftposp1should be negative.
Is the manual of the code wrong? Or my understanding is wrong somewhere?
the code is correct but the equation for the particle divergence (Eq. 22) in the manual (as well as the reference paper of Lee et al. 2008) is not.
The SPH divergence operator in general is D=Sum[(massj/rhopj)*(dr)*grad(Wij)] where dr is r(j)-r(i). In DualSPHysics however the particle distance is dr'=r(i)-r(j) creating a different operator D'=Sum[(massj/rhopj)*(dr')*grad(Wij)]. So to get the correct divergence, we need the minus: D=-D'.
Eq.22 in the manual changed the particle distance but did not add the "-". Apologies for my mistake, we will correct it in the next release.