diff --git a/ComputeGridGrid.cuh b/ComputeGridGrid.cuh
index 42428413c64651f7c2f401fd202224071ef31e31..12f9c2fae8d9f9a0f85a8dd4463bb2ba3277b913 100644
--- a/ComputeGridGrid.cuh
+++ b/ComputeGridGrid.cuh
@@ -8,7 +8,7 @@ __global__
void computeGridGridForce(const RigidBodyGrid* rho, const RigidBodyGrid* u,
Matrix3 basis_rho, Vector3 origin_rho,
Matrix3 basis_u, Vector3 origin_u,
- Vector3 * retForce, Vector3 * retTorque) {
+ Vector3 * retForce, Vector3 * retTorque, int gridNum) {
// printf("ComputeGridGridForce called\n");
/* extern __shared__ Vector3 force []; */
@@ -56,6 +56,11 @@ void computeGridGridForce(const RigidBodyGrid* rho, const RigidBodyGrid* u,
torque[tid] = t;
__syncthreads();
+ // debugging
+ float cutoff = 1e-3;
+ if (gridNum >= 0 && (abs(f.x) >= cutoff || abs(f.y) >= cutoff || abs(f.z) >= cutoff))
+ printf("GRIDFORCE: %d %f %f %f %f %f %f\n", gridNum, r_pos.x,r_pos.y,r_pos.z,f.x,f.y,f.z);
+
// Reduce force and torques
// http://www.cuvilib.com/Reduction.pdf
// RBTODO optimize
diff --git a/GrandBrownTown.cu b/GrandBrownTown.cu
index ca06f6a88af9afad839fb5c0ebc1e1df19cfba3b..92c6dfd694a7083938925a973a1271a977effdcb 100644
--- a/GrandBrownTown.cu
+++ b/GrandBrownTown.cu
@@ -399,7 +399,7 @@ void GrandBrownTown::run() {
int numBlocks = (num * numReplicas) / NUM_THREADS + (num * numReplicas % NUM_THREADS == 0 ? 0 : 1);
int tl = temperatureGrid.length();
- RBC.updateForces(); /* update RB forces before update particle positions... */
+ RBC.updateForces(s); /* update RB forces before update particle positions... */
// Call the kernel to update the positions of each particle
updateKernel<<< numBlocks, NUM_THREADS >>>(pos_d, forceInternal_d, type_d,
diff --git a/RigidBody.cu b/RigidBody.cu
index 1dd5f53fb6ee76f64d728b5999b8d619d80a342d..428859f2a1f4005f7c23fc5149d14a876d850f0a 100644
--- a/RigidBody.cu
+++ b/RigidBody.cu
@@ -77,7 +77,9 @@ void RigidBody::addLangevin(Vector3 w1, Vector3 w2) {
f = orientation * f; // return to lab frame
torq = orientation * torq;
-
+
+ // printf("LANGTORQUE: %f %f %f\n",torq.x,torq.y,torq.z);
+
addForce(f);
addTorque(torq);
}
diff --git a/RigidBody.h b/RigidBody.h
index 262fcf2ef328117cb8e1fa9254639049c4cb3f97..b9df591710f634e0964d64356a22a44fef616711 100644
--- a/RigidBody.h
+++ b/RigidBody.h
@@ -47,7 +47,7 @@ class RigidBody { // host side representation of rigid bodies
HOST DEVICE inline Vector3 getPosition() const { return position; }
HOST DEVICE inline Matrix3 getOrientation() const { return orientation; }
- HOST DEVICE inline Matrix3 getBasis() const { return orientation; }
+ // HOST DEVICE inline Matrix3 getBasis() const { return orientation; }
HOST DEVICE inline BigReal getMass() const { return t->mass; }
HOST DEVICE inline Vector3 getVelocity() const { return momentum/t->mass; }
HOST DEVICE inline Vector3 getAngularVelocity() const {
@@ -56,12 +56,13 @@ class RigidBody { // host side representation of rigid bodies
angularMomentum.z / t->inertia.z );
}
bool langevin;
+ Vector3 torque; // lab frame (except in integrate())
private:
String key;
/* static const SimParameters * simParams; */
Vector3 position;
- Matrix3 orientation;
+ Matrix3 orientation; /* rotation that brings RB coordinates into the lab frame */
Vector3 momentum;
Vector3 angularMomentum; // angular momentum along corresponding principal axes
@@ -81,7 +82,6 @@ private:
RigidBodyType* t; /* RBTODO: const? */
float timestep;
Vector3 force; // lab frame
- Vector3 torque; // lab frame (except in integrate())
bool isFirstStep;
diff --git a/RigidBodyController.cu b/RigidBodyController.cu
index 60f5fca60f097b81a244d3af39647f650063628a..31863513e1a1577f72c811b3aaf0ac906a228b32 100644
--- a/RigidBodyController.cu
+++ b/RigidBodyController.cu
@@ -150,7 +150,7 @@ void RigidBodyController::initializeForcePairs() {
}
}
-void RigidBodyController::updateForces() {
+void RigidBodyController::updateForces(int s) {
/*––{ RBTODO }–––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––.
| probably coalesce kernel calls, or move this to a device kernel caller |
| |
@@ -179,11 +179,25 @@ void RigidBodyController::updateForces() {
}
for (int i=0; i < forcePairs.size(); i++) {
- forcePairs[i].updateForces();
+ forcePairs[i].updateForces(i,s);
}
// RBTODO: see if there is a better way to sync
gpuErrchk(cudaDeviceSynchronize());
+
+ // debug
+ if (s %10 == 0) {
+ int tmp = 0;
+ for (int i = 0; i < rigidBodyByType.size(); i++) {
+ for (int j = 0; j < rigidBodyByType[i].size(); j++) {
+ RigidBody& rb = rigidBodyByType[i][j];
+ tmp++;
+ Vector3 p = rb.getPosition();
+ Vector3 t = rb.torque;
+ printf("RBTORQUE: %d %f %f %f %f %f %f\n", tmp, p.x, p.y, p.z, t.x,t.y,t.z);
+ }
+ }
+ }
}
void RigidBodyController::integrate(int step) {
// tell RBs to integrate
@@ -244,11 +258,14 @@ void RigidBodyController::integrate(int step) {
}
// RBTODO: bundle several rigidbodypair evaluations in single kernel call
-void RigidBodyForcePair::updateForces() {
+void RigidBodyForcePair::updateForces(int pairId, int s) {
// get the force/torque between a pair of rigid bodies
/* printf(" Updating rbPair forces\n"); */
const int numGrids = gridKeyId1.size();
+ if (s%10 != 0)
+ pairId = -1000;
+
// RBTODO: precompute certain common transformations and pass in kernel call
for (int i = 0; i < numGrids; i++) {
const int nb = numBlocks[i];
@@ -271,22 +288,22 @@ void RigidBodyForcePair::updateForces() {
| r = B'.ijk + c' |
`––––––––––––––––––./
*/
- Matrix3 B1 = rb1->getBasis()*type1->densityGrids[k1].getBasis();
- Vector3 c1 = rb1->getBasis()*type1->densityGrids[k1].getOrigin() + rb1->getPosition();
+ Matrix3 B1 = rb1->getOrientation()*type1->densityGrids[k1].getBasis();
+ Vector3 c1 = rb1->getOrientation()*type1->densityGrids[k1].getOrigin() + rb1->getPosition();
Matrix3 B2;
Vector3 c2;
/* printf(" Calculating grid forces\n"); */
if (!isPmf) { /* pair of RBs */
- B2 = rb2->getBasis()*type2->potentialGrids[k2].getBasis();
- c2 = rb2->getBasis()*type2->potentialGrids[k2].getOrigin() + rb2->getPosition();
+ B2 = rb2->getOrientation()*type2->potentialGrids[k2].getBasis();
+ c2 = rb2->getOrientation()*type2->potentialGrids[k2].getOrigin() + rb2->getPosition();
// RBTODO: get energy
computeGridGridForce<<< nb, numThreads >>>
(type1->rawDensityGrids_d[k1], type2->rawPotentialGrids_d[k2],
B1, c1, B2, c2,
- forces_d[i], torques_d[i]);
+ forces_d[i], torques_d[i], pairId+i);
} else { /* RB with a PMF */
// B2 = type2->rawPmfs[i].getBasis(); // not 100% certain k2 should be used rather than i
/// c2 = type2->rawPmfs[i].getOrigin();
@@ -296,7 +313,7 @@ void RigidBodyForcePair::updateForces() {
computeGridGridForce<<< nb, numThreads >>>
(type1->rawDensityGrids_d[k1], type2->rawPmfs_d[k2],
B1, c1, B2, c2,
- forces_d[i], torques_d[i]);
+ forces_d[i], torques_d[i], pairId+i);
}
}
@@ -326,11 +343,11 @@ void RigidBodyForcePair::updateForces() {
// transform torque from lab-frame origin to rb centers
// add forces to rbs
- Vector3 tmp;
- /* tmp = rb1->position; */
- /* printf("rb1->position: (%f,%f,%f)\n", tmp.x, tmp.y, tmp.z); */
- tmp = rb1->getPosition();
- printf("rb1->getPosition(): (%f,%f,%f)\n", tmp.x, tmp.y, tmp.z);
+ /* Vector3 tmp; */
+ /* /\* tmp = rb1->position; *\/ */
+ /* /\* printf("rb1->position: (%f,%f,%f)\n", tmp.x, tmp.y, tmp.z); *\/ */
+ /* tmp = rb1->getPosition(); */
+ /* printf("rb1->getPosition(): (%f,%f,%f)\n", tmp.x, tmp.y, tmp.z); */
Vector3 t1 = t - rb1->getPosition().cross( f );
rb1->addForce( f);
rb1->addTorque(t1);
@@ -509,6 +526,7 @@ void RigidBodyController::print(int step) {
trajFile << "# RigidBody trajectory file" << std::endl;
printLegend(trajFile);
}
+ printf("WRITING RIGID BODY COORDINATES AT STEP %d\n",step);
printData(step,trajFile);
trajFile.flush();
}
@@ -596,7 +614,6 @@ void RigidBodyController::printLegend(std::ofstream &file) {
<< " angVelX angVelY angVelZ" << std::endl;
}
void RigidBodyController::printData(int step,std::ofstream &file) {
- printf("WRITING RIGID BODY COORDINATES AT STEP %d\n",step);
// tell RBs to integrate
for (int i = 0; i < rigidBodyByType.size(); i++) {
for (int j = 0; j < rigidBodyByType[i].size(); j++) {
diff --git a/RigidBodyController.h b/RigidBodyController.h
index 6e6861493284eef5ad55fc1e9d6a081e2a57f7e6..4acfa5f4ff1b8ff1f6da4f6ee1714e17733a7cc6 100644
--- a/RigidBodyController.h
+++ b/RigidBodyController.h
@@ -60,7 +60,7 @@ private:
std::vector<Vector3*> torques;
std::vector<Vector3*> torques_d;
- void updateForces();
+ void updateForces(int pairId, int s);
};
class RigidBodyController {
@@ -71,7 +71,7 @@ public:
RigidBodyController(const Configuration& c, const char* outArg);
void integrate(int step);
- void updateForces();
+ void updateForces(int s);
private:
void copyGridsToDevice();