diff --git a/ComputeGridGrid.cuh b/ComputeGridGrid.cuh
index 55e2bc4f7ef45f91dba720f12eb861ed5bb4871e..8befee02e254edf16e4760c8e05e3afca80be55b 100644
--- a/ComputeGridGrid.cuh
+++ b/ComputeGridGrid.cuh
@@ -68,8 +68,8 @@ void computeGridGridForce(const RigidBodyGrid* rho, const RigidBodyGrid* u,
if (tid == 0) {
retForce[blockIdx.x] = force[0];
retTorque[blockIdx.x] = torque[0];
- printf("GPU force: (%f,%f,%f)\n", force[0].x, force[0].y, force[0].z);
- printf("GPU force0: (%f,%f,%f)\n", f.x, f.y, f.z);
+ /* printf("GPU force: (%f,%f,%f)\n", force[0].x, force[0].y, force[0].z); */
+ /* printf("GPU force0: (%f,%f,%f)\n", f.x, f.y, f.z); */
}
}
diff --git a/Configuration.cpp b/Configuration.cpp
index 3a759924bf421f0f956e0d9f715b2ed314894fed..258ce8170aadd486d702fe7ff01ac721cb0fcb97 100644
--- a/Configuration.cpp
+++ b/Configuration.cpp
@@ -664,8 +664,6 @@ int Configuration::readParameters(const char * config_file) {
part[++currPart] = BrownianParticleType(value);
currPartClass = partClassPart;
}
- else if (param == String("gridFile"))
- partGridFile[currPart] = value;
else if (param == String("forceXGridFile"))
partForceXGridFile[currPart] = value;
else if (param == String("forceYGridFile"))
@@ -784,9 +782,15 @@ int Configuration::readParameters(const char * config_file) {
else if (param == String("num")) {
if (currPartClass == partClassPart)
part[currPart].num = atoi(value.val());
- else if (currPartClass == partClassRB)
+ else if (currPartClass == partClassRB)
rigidBody[currRB].num = atoi(value.val());
}
+ else if (param == String("gridFile")) {
+ if (currPartClass == partClassPart)
+ partGridFile[currPart] = value;
+ else if (currPartClass == partClassRB)
+ rigidBody[currRB].addPMF(value);
+ }
// UNKNOWN
else {
printf("WARNING: Unrecognized keyword `%s'.\n", param.val());
diff --git a/RigidBody.cu b/RigidBody.cu
index 9a4b1a406272027fee682c63c792425477e423cf..2a6347876e12c5b7e735f6773ab6f940b8740e4d 100644
--- a/RigidBody.cu
+++ b/RigidBody.cu
@@ -92,7 +92,7 @@ void RigidBody::integrate(int startFinishAll) {
Vector3 trans; // = *p_trans;
Matrix3 rot = Matrix3(1); // = *p_rot;
- printf("Rigid Body force\n");
+ /* printf("Rigid Body force\n"); */
#ifdef DEBUGM
switch (startFinishAll) {
diff --git a/RigidBodyController.cu b/RigidBodyController.cu
index 253ec54e81975358f2c31efadd94a92571750a06..7f11643bf6a193c477b5868efc67eb569285ba47 100644
--- a/RigidBodyController.cu
+++ b/RigidBodyController.cu
@@ -105,7 +105,7 @@ void RigidBodyController::initializeForcePairs() {
RigidBody* rb2 = &(rbs2[j]);
printf(" pushing RB force pair for %d:%d\n",i,j);
- RigidBodyForcePair fp = RigidBodyForcePair(&(t1),&(t2),rb1,rb2,gridKeyId1,gridKeyId2);
+ RigidBodyForcePair fp = RigidBodyForcePair(&(t1),&(t2),rb1,rb2,gridKeyId1,gridKeyId2, false);
gpuErrchk(cudaDeviceSynchronize()); /* RBTODO: this should be extraneous */
forcePairs.push_back( fp );
printf(" done pushing RB force pair for %d:%d\n",i,j);
@@ -114,6 +114,40 @@ void RigidBodyController::initializeForcePairs() {
}
}
}
+
+ // add Pmfs (not a true pairwise RB interaction; hacky implementation)
+ for (int ti = 0; ti < conf.numRigidTypes; ti++) {
+ RigidBodyType& t1 = conf.rigidBody[ti];
+
+ const std::vector<String>& keys1 = t1.densityGridKeys;
+ const std::vector<String>& keys2 = t1.pmfKeys;
+ std::vector<int> gridKeyId1;
+ std::vector<int> gridKeyId2;
+
+ // Loop over all pairs of grid keys (e.g. "Elec")
+ bool paired = false;
+ for(int k1 = 0; k1 < keys1.size(); k1++) {
+ for(int k2 = 0; k2 < keys2.size(); k2++) {
+ if ( keys1[k1] == keys2[k2] ) {
+ gridKeyId1.push_back(k1);
+ gridKeyId2.push_back(k2);
+ paired = true;
+ }
+ }
+ }
+ if (paired) {
+ // found matching keys => calculate force between all grid pairs
+ std::vector<RigidBody>& rbs1 = rigidBodyByType[ti];
+
+ // Loop over rigid bodies of these types
+ for (int i = 0; i < rbs1.size(); i++) {
+ RigidBody* rb1 = &(rbs1[i]);
+ RigidBodyForcePair fp = RigidBodyForcePair(&(t1),&(t1),rb1,rb1,gridKeyId1,gridKeyId2, true);
+ gpuErrchk(cudaDeviceSynchronize()); /* RBTODO: this should be extraneous */
+ forcePairs.push_back( fp );
+ }
+ }
+ }
}
void RigidBodyController::updateForces() {
@@ -138,11 +172,9 @@ void RigidBodyController::updateForces() {
for (int i=0; i < forcePairs.size(); i++) {
forcePairs[i].updateForces();
}
- // get 3rd law forces and torques
-
+
// RBTODO: see if there is a better way to sync
gpuErrchk(cudaDeviceSynchronize());
-
}
void RigidBodyController::integrate(int step) {
// tell RBs to integrate
@@ -205,7 +237,7 @@ void RigidBodyController::integrate(int step) {
// RBTODO: bundle several rigidbodypair evaluations in single kernel call
void RigidBodyForcePair::updateForces() {
// get the force/torque between a pair of rigid bodies
- printf(" Updating rbPair forces\n");
+ /* printf(" Updating rbPair forces\n"); */
const int numGrids = gridKeyId1.size();
// RBTODO: precompute certain common transformations and pass in kernel call
@@ -214,16 +246,25 @@ void RigidBodyForcePair::updateForces() {
const int k1 = gridKeyId1[i];
const int k2 = gridKeyId2[i];
- printf(" Calculating grid forces\n");
-
- // RBTODO: add energy
- computeGridGridForce<<< nb, numThreads >>>
- (type1->rawDensityGrids_d[k1], type2->rawPotentialGrids_d[k2],
- rb1->getBasis(), rb1->getPosition(), /* RBTODO: include offset from grid */
- rb2->getBasis(), rb2->getPosition(),
- forces_d[i], torques_d[i]);
-
+ /* printf(" Calculating grid forces\n"); */
+ if (!isPmf) { /* pair of RBs */
+ // RBTODO: add energy
+ computeGridGridForce<<< nb, numThreads >>>
+ (type1->rawDensityGrids_d[k1], type2->rawPotentialGrids_d[k2],
+ rb1->getBasis(), rb1->getPosition(), /* RBTODO: include offset from grid */
+ rb2->getBasis(), rb2->getPosition(),
+ forces_d[i], torques_d[i]);
+ } else { /* RB with a PMF */
+ computeGridGridForce<<< nb, numThreads >>>
+ (type1->rawDensityGrids_d[k1], type2->rawPmfs_d[k2],
+ rb1->getBasis(), rb1->getPosition(), /* RBTODO: include offset from grid */
+ type2->rawPmfs[i].getBasis(), type2->rawPmfs[i].getOrigin(),
+ forces_d[i], torques_d[i]);
+ }
+
}
+
+ // RBTODO better way to sync?
gpuErrchk(cudaDeviceSynchronize());
for (int i = 0; i < numGrids; i++) {
const int nb = numBlocks[i];
@@ -232,7 +273,6 @@ void RigidBodyForcePair::updateForces() {
gpuErrchk(cudaMemcpy(torques[i], torques_d[i], sizeof(Vector3)*nb,
cudaMemcpyDeviceToHost));
}
-
gpuErrchk(cudaDeviceSynchronize());
// sum forces + torques
@@ -248,17 +288,19 @@ void RigidBodyForcePair::updateForces() {
}
// transform torque from lab-frame origin to rb centers
- Vector3 t1 = t - rb1->getPosition().cross( f );
- Vector3 t2 = -t - rb2->getPosition().cross( -f );
-
// add forces to rbs
+ Vector3 t1 = t - rb1->getPosition().cross( f );
rb1->addForce( f);
rb1->addTorque(t1);
- rb2->addForce(-f);
- rb2->addTorque(t2);
- printf("force: (%f,%f,%f)\n",f.x,f.y,f.z);
- printf("torque: (%f,%f,%f)\n",t1.x,t1.y,t1.z);
+ if (! isPmf) {
+ Vector3 t2 = -t - rb2->getPosition().cross( -f );
+ rb2->addForce(-f);
+ rb2->addTorque(t2);
+ }
+
+ /* printf("force: (%f,%f,%f)\n",f.x,f.y,f.z); */
+ /* printf("torque: (%f,%f,%f)\n",t1.x,t1.y,t1.z); */
}
void RigidBodyController::copyGridsToDevice() {
@@ -350,15 +392,49 @@ void RigidBodyController::copyGridsToDevice() {
}
}
+ for (int i = 0; i < conf.numRigidTypes; i++) {
+ printf("Copying PMFs for RB %d\n",i);
+ RigidBodyType& rb = conf.rigidBody[i];
+
+ int ng = rb.numPmfs;
+ rb.rawPmfs_d = new RigidBodyGrid*[ng]; /* not 100% sure this is needed, possible memory leak */
+
+ printf(" RigidBodyType %d: numPmfs = %d\n", i, ng);
+
+ // copy pmf grid data to device
+ for (int gid = 0; gid < ng; gid++) {
+ RigidBodyGrid g = rb.rawPmfs[gid];
+ int len = g.getSize();
+ float* tmpData;
+ // tmpData = new float*[len];
+
+ size_t sz = sizeof(RigidBodyGrid);
+ gpuErrchk(cudaMalloc((void **) &(rb.rawPmfs_d[gid]), sz));
+ gpuErrchk(cudaMemcpy( rb.rawPmfs_d[gid], &g,
+ sz, cudaMemcpyHostToDevice ));
+
+ // allocate grid data on device
+ // copy temporary host pointer to device pointer
+ // copy data to device through temporary host pointer
+ sz = sizeof(float) * len;
+ gpuErrchk(cudaMalloc((void **) &tmpData, sz));
+ // sz = sizeof(float) * len;
+ gpuErrchk(cudaMemcpy( tmpData, rb.rawPmfs[gid].val, sz, cudaMemcpyHostToDevice));
+ gpuErrchk(cudaMemcpy( &(rb.rawPmfs_d[gid]->val), &tmpData,
+ sizeof(float*), cudaMemcpyHostToDevice));
+
+ }
+ }
+
gpuErrchk(cudaDeviceSynchronize());
printf("Done copying RBs\n");
- // DEBUG
- RigidBodyType& rb = conf.rigidBody[0];
- printRigidBodyGrid<<<1,1>>>( rb.rawPotentialGrids_d[0] );
- gpuErrchk(cudaDeviceSynchronize());
- printRigidBodyGrid<<<1,1>>>( rb.rawDensityGrids_d[0] );
- gpuErrchk(cudaDeviceSynchronize());
+ /* // DEBUG */
+ /* RigidBodyType& rb = conf.rigidBody[0]; */
+ /* printRigidBodyGrid<<<1,1>>>( rb.rawPotentialGrids_d[0] ); */
+ /* gpuErrchk(cudaDeviceSynchronize()); */
+ /* printRigidBodyGrid<<<1,1>>>( rb.rawDensityGrids_d[0] ); */
+ /* gpuErrchk(cudaDeviceSynchronize()); */
}
void RigidBodyController::print(int step) {
diff --git a/RigidBodyController.h b/RigidBodyController.h
index 5bfde4eea9ceef057727eb8a0119e21067d41ba5..6e6861493284eef5ad55fc1e9d6a081e2a57f7e6 100644
--- a/RigidBodyController.h
+++ b/RigidBodyController.h
@@ -16,9 +16,9 @@ class RigidBodyForcePair {
public:
RigidBodyForcePair(RigidBodyType* t1, RigidBodyType* t2,
RigidBody* rb1, RigidBody* rb2,
- std::vector<int> gridKeyId1, std::vector<int> gridKeyId2) :
+ std::vector<int> gridKeyId1, std::vector<int> gridKeyId2, bool isPmf) :
type1(t1), type2(t2), rb1(rb1), rb2(rb2),
- gridKeyId1(gridKeyId1), gridKeyId2(gridKeyId2)
+ gridKeyId1(gridKeyId1), gridKeyId2(gridKeyId2), isPmf(isPmf)
{
printf(" Constructing RB force pair...\n");
initialize();
@@ -26,7 +26,7 @@ public:
}
RigidBodyForcePair(const RigidBodyForcePair& o) :
type1(o.type1), type2(o.type2), rb1(o.rb1), rb2(o.rb2),
- gridKeyId1(o.gridKeyId1), gridKeyId2(o.gridKeyId2) {
+ gridKeyId1(o.gridKeyId1), gridKeyId2(o.gridKeyId2), isPmf(o.isPmf) {
printf(" Copying RB force pair...\n");
initialize();
}
@@ -44,6 +44,8 @@ private:
static const int numThreads = NUMTHREADS;
+ bool isPmf;
+
RigidBodyType* type1;
RigidBodyType* type2;
RigidBody* rb1;
diff --git a/RigidBodyType.cu b/RigidBodyType.cu
index 2d9ae682e99025ba44e8a8322e96a07803b7a633..bf0f416793c8f7b9c9e1aa33eafd37220a26607e 100644
--- a/RigidBodyType.cu
+++ b/RigidBodyType.cu
@@ -10,10 +10,12 @@ void RigidBodyType::clear() {
// TODO: make sure that this actually removes grid data
potentialGrids.clear();
densityGrids.clear();
-
+ pmfs.clear();
+
potentialGridKeys.clear();
densityGridKeys.clear();
-
+ pmfKeys.clear();
+
if (numPotGrids > 0) delete[] rawPotentialGrids;
if (numDenGrids > 0) delete[] rawDensityGrids;
rawPotentialGrids = NULL;
@@ -126,12 +128,30 @@ void RigidBodyType::addDensityGrid(String s) {
densityGrids.push_back( g );
densityGridKeys.push_back( key );
}
+void RigidBodyType::addPMF(String s) {
+ // tokenize and return
+ int numTokens = s.tokenCount();
+ if (numTokens != 2) {
+ printf("ERROR: could not add Grid.\n"); // TODO improve this message
+ exit(1);
+ }
+ String* token = new String[numTokens];
+ s.tokenize(token);
+ String key = token[0];
+ BaseGrid g(token[1]);
+
+ pmfs.push_back( g );
+ pmfKeys.push_back( key );
+}
void RigidBodyType::updateRaw() {
if (numPotGrids > 0) delete[] rawPotentialGrids;
if (numDenGrids > 0) delete[] rawDensityGrids;
+ if (numDenGrids > 0) delete[] rawPmfs;
numPotGrids = potentialGrids.size();
numDenGrids = densityGrids.size();
+ numPmfs = pmfs.size();
+
if (numPotGrids > 0) {
rawPotentialGrids = new RigidBodyGrid[numPotGrids];
rawPotentialBases = new Matrix3[numPotGrids];
@@ -142,6 +162,8 @@ void RigidBodyType::updateRaw() {
rawDensityBases = new Matrix3[numDenGrids];
rawDensityOrigins = new Vector3[numDenGrids];
}
+ if (numPmfs > 0)
+ rawPmfs = new BaseGrid[numPmfs];
for (int i=0; i < numPotGrids; i++) {
rawPotentialGrids[i] = potentialGrids[i];
@@ -153,6 +175,8 @@ void RigidBodyType::updateRaw() {
rawDensityBases[i] = densityGrids[i].getBasis();
rawDensityOrigins[i] = densityGrids[i].getOrigin();
}
+ for (int i=0; i < numPmfs; i++)
+ rawPmfs[i] = pmfs[i];
}
diff --git a/RigidBodyType.h b/RigidBodyType.h
index 707b3d1c60b5c9301c30523ecd12e4e571285d2a..c9e32586aa84b09b9c45c8ebaa2544182df0d708 100644
--- a/RigidBodyType.h
+++ b/RigidBodyType.h
@@ -28,7 +28,7 @@ public:
RigidBodyType(const String& name = "") :
name(name), num(0),
reservoir(NULL), mass(1.0f), inertia(), transDamping(),
- rotDamping(), numPotGrids(0), numDenGrids(0) { }
+ rotDamping(), numPotGrids(0), numDenGrids(0), numPmfs(0) { }
/* RigidBodyType(const RigidBodyType& src) { copy(src); } */
~RigidBodyType() { clear(); }
@@ -37,6 +37,7 @@ RigidBodyType(const String& name = "") :
void addPotentialGrid(String s);
void addDensityGrid(String s);
+ void addPMF(String s);
void updateRaw();
void setDampingCoeffs(float timestep);
@@ -57,16 +58,20 @@ public:
std::vector<String> potentialGridKeys;
std::vector<String> densityGridKeys;
+ std::vector<String> pmfKeys;
std::vector<BaseGrid> potentialGrids;
std::vector<BaseGrid> densityGrids;
+ std::vector<BaseGrid> pmfs;
// RBTODO: clear std::vectors after initialization
// duplicates of std::vector grids for device
int numPotGrids;
int numDenGrids;
+ int numPmfs;
RigidBodyGrid* rawPotentialGrids;
RigidBodyGrid* rawDensityGrids;
+ BaseGrid* rawPmfs;
Matrix3* rawPotentialBases;
Matrix3* rawDensityBases;
Vector3* rawPotentialOrigins;
@@ -75,6 +80,7 @@ public:
// device pointers
RigidBodyGrid** rawPotentialGrids_d;
RigidBodyGrid** rawDensityGrids_d;
+ RigidBodyGrid** rawPmfs_d;
};
diff --git a/notes.org b/notes.org
index b64defb2f7403d3592df54c34116f451c54f9b0c..a9e9349214cf5c89a31d8a169dd4e1214c0ade9e 100644
--- a/notes.org
+++ b/notes.org
@@ -3,8 +3,10 @@
** branch rb-device: have RBType objects hold device pointers for raw RigidBodyGrids
* TODO eventually
-** read restart file
** multiple "PMF" grids associated with each RB type
+** read restart file
+** throw exception if an RB type uses a grid key multiple times
+** improve efficiency of force evaluation kernel
** each RB gets temperature from grid
** optionally don't eval forces every ts