Skip to content
Snippets Groups Projects
Commit ff11ff2b authored by cmaffeo2's avatar cmaffeo2
Browse files

added pmfs for rigid bodies!

parent f891f45d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
ComputeGridGrid.cuh
+ 2
− 2
View file @ ff11ff2b
......@@ -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); */
}
}
......
Configuration.cpp
+ 7
− 3
View file @ ff11ff2b
......@@ -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());
......
RigidBody.cu
+ 1
− 1
View file @ ff11ff2b
......@@ -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) {
......
RigidBodyController.cu
+ 104
− 28
View file @ ff11ff2b
......@@ -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) {
......
RigidBodyController.h
+ 5
− 3
View file @ ff11ff2b
......@@ -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;
......
RigidBodyType.cu
+ 26
− 2
View file @ ff11ff2b
......@@ -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];
}
RigidBodyType.h
+ 7
− 1
View file @ ff11ff2b
......@@ -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;
};
notes.org
+ 3
− 1
View file @ ff11ff2b
......@@ -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
......
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