Skip to content

GitLab

Commit fd792ec1 authored by Robert T. McGibbon's avatar Robert T. McGibbon
Browse files

Fix for r=0 in periodicdistance

parent c741523e
Hide whitespace changes
Inline Side-by-side
Showing with 18 additions and 14 deletions
platforms/cuda/src/CudaExpressionUtilities.cpp
View file @ fd792ec1
...@@ -109,7 +109,8 @@ void CudaExpressionUtilities::processExpression(stringstream& out, const Express ...@@ -109,7 +109,8 @@ void CudaExpressionUtilities::processExpression(stringstream& out, const Express
} }
out << ");\n"; out << ");\n";
out << "APPLY_PERIODIC_TO_DELTA(periodicDistance_delta)\n"; out << "APPLY_PERIODIC_TO_DELTA(periodicDistance_delta)\n";
out << tempType << " periodicDistance_rinv = RSQRT(periodicDistance_delta.x*periodicDistance_delta.x + periodicDistance_delta.y*periodicDistance_delta.y + periodicDistance_delta.z*periodicDistance_delta.z);\n"; out << tempType << " periodicDistance_r2 = periodicDistance_delta.x*periodicDistance_delta.x + periodicDistance_delta.y*periodicDistance_delta.y + periodicDistance_delta.z*periodicDistance_delta.z;\n";
out << tempType << " periodicDistance_rinv = RSQRT(periodicDistance_r2);\n";
for (int j = 0; j < nodes.size(); j++) { for (int j = 0; j < nodes.size(); j++) {
const vector<int>& derivOrder = dynamic_cast<const Operation::Custom*>(&nodes[j]->getOperation())->getDerivOrder(); const vector<int>& derivOrder = dynamic_cast<const Operation::Custom*>(&nodes[j]->getOperation())->getDerivOrder();
int argIndex = -1; int argIndex = -1;
...@@ -123,17 +124,17 @@ void CudaExpressionUtilities::processExpression(stringstream& out, const Express ...@@ -123,17 +124,17 @@ void CudaExpressionUtilities::processExpression(stringstream& out, const Express
if (argIndex == -1) if (argIndex == -1)
out << nodeNames[j] << " = RECIP(periodicDistance_rinv);\n"; out << nodeNames[j] << " = RECIP(periodicDistance_rinv);\n";
else if (argIndex == 0) else if (argIndex == 0)
out << nodeNames[j] << " = periodicDistance_delta.x*periodicDistance_rinv;\n"; out << nodeNames[j] << " = periodicDistance_delta.x*periodicDistance_rinv*(periodicDistance_r2>0);\n";
else if (argIndex == 1) else if (argIndex == 1)
out << nodeNames[j] << " = periodicDistance_delta.y*periodicDistance_rinv;\n"; out << nodeNames[j] << " = periodicDistance_delta.y*periodicDistance_rinv*(periodicDistance_r2>0);\n";
else if (argIndex == 2) else if (argIndex == 2)
out << nodeNames[j] << " = periodicDistance_delta.z*periodicDistance_rinv;\n"; out << nodeNames[j] << " = periodicDistance_delta.z*periodicDistance_rinv*(periodicDistance_r2>0);\n";
else if (argIndex == 3) else if (argIndex == 3)
out << nodeNames[j] << " = -periodicDistance_delta.x*periodicDistance_rinv;\n"; out << nodeNames[j] << " = -periodicDistance_delta.x*periodicDistance_rinv*(periodicDistance_r2>0);\n";
else if (argIndex == 4) else if (argIndex == 4)
out << nodeNames[j] << " = -periodicDistance_delta.y*periodicDistance_rinv;\n"; out << nodeNames[j] << " = -periodicDistance_delta.y*periodicDistance_rinv*(periodicDistance_r2>0);\n";
else if (argIndex == 5) else if (argIndex == 5)
out << nodeNames[j] << " = -periodicDistance_delta.z*periodicDistance_rinv;\n"; out << nodeNames[j] << " = -periodicDistance_delta.z*periodicDistance_rinv*(periodicDistance_r2>0);\n";
} }
} }
else { else {
......
platforms/opencl/src/OpenCLExpressionUtilities.cpp
View file @ fd792ec1
...@@ -109,7 +109,8 @@ void OpenCLExpressionUtilities::processExpression(stringstream& out, const Expre ...@@ -109,7 +109,8 @@ void OpenCLExpressionUtilities::processExpression(stringstream& out, const Expre
} }
out << ");\n"; out << ");\n";
out << "APPLY_PERIODIC_TO_DELTA(periodicDistance_delta)\n"; out << "APPLY_PERIODIC_TO_DELTA(periodicDistance_delta)\n";
out << tempType << " periodicDistance_rinv = RSQRT(periodicDistance_delta.x*periodicDistance_delta.x + periodicDistance_delta.y*periodicDistance_delta.y + periodicDistance_delta.z*periodicDistance_delta.z);\n"; out << tempType << " periodicDistance_r2 = periodicDistance_delta.x*periodicDistance_delta.x + periodicDistance_delta.y*periodicDistance_delta.y + periodicDistance_delta.z*periodicDistance_delta.z;\n";
out << tempType << " periodicDistance_rinv = RSQRT(periodicDistance_r2);\n";
for (int j = 0; j < nodes.size(); j++) { for (int j = 0; j < nodes.size(); j++) {
const vector<int>& derivOrder = dynamic_cast<const Operation::Custom*>(&nodes[j]->getOperation())->getDerivOrder(); const vector<int>& derivOrder = dynamic_cast<const Operation::Custom*>(&nodes[j]->getOperation())->getDerivOrder();
int argIndex = -1; int argIndex = -1;
...@@ -123,17 +124,17 @@ void OpenCLExpressionUtilities::processExpression(stringstream& out, const Expre ...@@ -123,17 +124,17 @@ void OpenCLExpressionUtilities::processExpression(stringstream& out, const Expre
if (argIndex == -1) if (argIndex == -1)
out << nodeNames[j] << " = RECIP(periodicDistance_rinv);\n"; out << nodeNames[j] << " = RECIP(periodicDistance_rinv);\n";
else if (argIndex == 0) else if (argIndex == 0)
out << nodeNames[j] << " = periodicDistance_delta.x*periodicDistance_rinv;\n"; out << nodeNames[j] << " = periodicDistance_delta.x*periodicDistance_rinv*(periodicDistance_r2>0);\n";
else if (argIndex == 1) else if (argIndex == 1)
out << nodeNames[j] << " = periodicDistance_delta.y*periodicDistance_rinv;\n"; out << nodeNames[j] << " = periodicDistance_delta.y*periodicDistance_rinv*(periodicDistance_r2>0);\n";
else if (argIndex == 2) else if (argIndex == 2)
out << nodeNames[j] << " = periodicDistance_delta.z*periodicDistance_rinv;\n"; out << nodeNames[j] << " = periodicDistance_delta.z*periodicDistance_rinv*(periodicDistance_r2>0);\n";
else if (argIndex == 3) else if (argIndex == 3)
out << nodeNames[j] << " = -periodicDistance_delta.x*periodicDistance_rinv;\n"; out << nodeNames[j] << " = -periodicDistance_delta.x*periodicDistance_rinv*(periodicDistance_r2>0);\n";
else if (argIndex == 4) else if (argIndex == 4)
out << nodeNames[j] << " = -periodicDistance_delta.y*periodicDistance_rinv;\n"; out << nodeNames[j] << " = -periodicDistance_delta.y*periodicDistance_rinv*(periodicDistance_r2>0);\n";
else if (argIndex == 5) else if (argIndex == 5)
out << nodeNames[j] << " = -periodicDistance_delta.z*periodicDistance_rinv;\n"; out << nodeNames[j] << " = -periodicDistance_delta.z*periodicDistance_rinv*(periodicDistance_r2>0);\n";
} }
} }
else { else {
......
platforms/reference/src/ReferenceKernels.cpp
View file @ fd792ec1
...@@ -1488,6 +1488,8 @@ double ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::evaluat ...@@ -1488,6 +1488,8 @@ double ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::evaluat
delta -= boxVectors[1]*floor(delta[1]/boxVectors[1][1]+0.5); delta -= boxVectors[1]*floor(delta[1]/boxVectors[1][1]+0.5);
delta -= boxVectors[0]*floor(delta[0]/boxVectors[0][0]+0.5); delta -= boxVectors[0]*floor(delta[0]/boxVectors[0][0]+0.5);
double r = sqrt(delta.dot(delta)); double r = sqrt(delta.dot(delta));
if (r == 0)
return 0.0;
if (argIndex < 3) if (argIndex < 3)
return delta[argIndex]/r; return delta[argIndex]/r;
return -delta[argIndex-3]/r; return -delta[argIndex-3]/r;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment