diff --git a/Configuration.cpp b/Configuration.cpp
index 2efe3e77603d7e8a54422e7ce7609efb94111578..42b9f27ea0e49fb4c98cde409c47d4aebfcb2684 100644
--- a/Configuration.cpp
+++ b/Configuration.cpp
@@ -295,6 +295,15 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
}
}
+ // Count number of particles of each type
+ numPartsOfType = new int[numParts];
+ for (int i = 0; i < numParts; ++i) {
+ numPartsOfType[i] = 0;
+ }
+ for (int i = 0; i < num; ++i) {
+ ++numPartsOfType[type[i]];
+ }
+
// Some geometric stuff that should be gotten rid of.
Vector3 buffer = (sys->getCenter() + 2.0f*sys->getOrigin())/3.0f;
initialZ = buffer.z;
@@ -381,6 +390,7 @@ Configuration::~Configuration() {
delete[] partForceZGridFile;
delete[] partDiffusionGridFile;
delete[] partReservoirFile;
+ delete[] partRigidBodyGrid;
// TODO: plug memory leaks
if (partsFromFile != NULL) delete[] partsFromFile;
@@ -391,6 +401,8 @@ Configuration::~Configuration() {
if (angles != NULL) delete[] angles;
if (dihedrals != NULL) delete[] dihedrals;
+ delete[] numPartsOfType;
+
// Table parameters
delete[] partTableFile;
delete[] partTableIndex0;
@@ -588,13 +600,14 @@ int Configuration::readParameters(const char * config_file) {
partForceZGridFile = new String[numParts];
partDiffusionGridFile = new String[numParts];
partReservoirFile = new String[numParts];
-
+ partRigidBodyGrid = new std::vector<String>[numParts];
+
// Allocate the table variables.
partTableFile = new String[numParts*numParts];
partTableIndex0 = new int[numParts*numParts];
partTableIndex1 = new int[numParts*numParts];
- // Allocate rigid body types
+ // Allocate rigid body types
rigidBody = new RigidBodyType[numRigidTypes];
int btfcap = 10;
@@ -764,7 +777,9 @@ int Configuration::readParameters(const char * config_file) {
}
if (readDihedralFile(value, ++currDihedral))
numTabDihedralFiles++;
- }
+ } else if (param == String("rigidBodyPotential")) {
+ partRigidBodyGrid[currPart].push_back(value);
+ }
// RIGID BODY
else if (param == String("rigidBody")) {
// part[++currPart] = BrownianParticleType(value);
diff --git a/Configuration.h b/Configuration.h
index 058b0bb31ecb2f4e39a472b42517a313d5b33c1a..9d23bb83bf6232c1b4a044ff77004aa155b2372e 100644
--- a/Configuration.h
+++ b/Configuration.h
@@ -156,6 +156,7 @@ public:
int numExcludes;
int numAngles;
int numDihedrals;
+ int* numPartsOfType;
String partFile;
String bondFile;
String excludeFile;
@@ -167,6 +168,7 @@ public:
bool readAnglesFromFile;
bool readDihedralsFromFile;
String* partGridFile;
+ std::vector<String>* partRigidBodyGrid;
String* partDiffusionGridFile;
String* partForceXGridFile;
String* partForceYGridFile;
diff --git a/GrandBrownTown.cu b/GrandBrownTown.cu
index 4eb25721ed2c5549c4b25d3fb8f03b4afd62f579..bc4788cd471cd789ca3e5b420c3152013b89635f 100644
--- a/GrandBrownTown.cu
+++ b/GrandBrownTown.cu
@@ -451,7 +451,7 @@ void GrandBrownTown::run() {
int numBlocks = (num * numReplicas) / NUM_THREADS + (num * numReplicas % NUM_THREADS == 0 ? 0 : 1);
int tl = temperatureGrid.length();
- RBC.updateForces(s); /* update RB forces before update particle positions... */
+ RBC.updateForces(internal->getPos_d(), internal->getForceInternal_d(), s); /* update RB forces before update particle positions... */
/* DEBUG: reduce net force on particles
{
diff --git a/RigidBodyController.cu b/RigidBodyController.cu
index e3744a2db8583ac62e8b59064db37f5ae3df110e..e7bc0c55477d9d59ad89efc69e6261a6502bdb99 100644
--- a/RigidBodyController.cu
+++ b/RigidBodyController.cu
@@ -59,7 +59,8 @@ RigidBodyController::RigidBodyController(const Configuration& c, const char* out
random = new RandomCPU(conf.seed + 1); /* +1 to avoid using same seed as RandomCUDA */
-initializeForcePairs();
+ initializeForcePairs();
+ initializeParticleLists();
}
RigidBodyController::~RigidBodyController() {
for (int i = 0; i < rigidBodyByType.size(); i++)
@@ -166,8 +167,90 @@ void RigidBodyController::initializeForcePairs() {
forcePairs[i].initialize();
}
+
+void RigidBodyController::initializeParticleLists() {
+ // Populate RigidBodyType.particles
+
+ // TODO: ensure no duplicates in conf.partRigidBodyGrid[i]
-void RigidBodyController::updateForces(int s) {
+ // Allocate RB type's numParticles array
+ for (int rb = 0; rb < conf.numRigidTypes; ++rb) {
+ RigidBodyType& t = conf.rigidBody[rb];
+ t.numParticles = new int[t.numPotGrids];
+ for (int i = 0; i < t.numPotGrids; ++i) t.numParticles[i] = 0;
+ }
+
+ // Count the number of particles; Loop over particle types
+ for (int i = 0; i < conf.numParts; ++i) {
+
+ // Loop over rigid body grid names associated with particle type
+ const std::vector<String>& gridNames = conf.partRigidBodyGrid[i];
+ for (int j = 0; j < gridNames.size(); ++j) {
+
+ // Loop over RB types
+ for (int rb = 0; rb < conf.numRigidTypes; ++rb) {
+ RigidBodyType& t = conf.rigidBody[rb];
+ const std::vector<String>& keys = t.potentialGridKeys;
+
+ // Loop over potential grids
+ for(int k = 0; k < keys.size(); k++) {
+ // printf(" checking grid keys ");
+ if (gridNames[j] == keys[k])
+ t.numParticles[k] += conf.numPartsOfType[i];
+ }
+ }
+ }
+ }
+
+ // Allocate each particles array
+ for (int rb = 0; rb < conf.numRigidTypes; ++rb) {
+ RigidBodyType& t = conf.rigidBody[rb];
+ t.particles = new int*[t.numPotGrids];
+ for (int i = 0; i < t.numPotGrids; ++i) {
+ t.particles[i] = new int[t.numParticles[i]];
+ t.numParticles[i] = 0; // now use this as a counter
+ }
+ }
+
+ // Set the number of particles; Loop over particle types
+ for (int i = 0; i < conf.numParts; ++i) {
+ int tmp[conf.numPartsOfType[i]]; // temporary array holding particles of type i
+ int currId = 0;
+ for (int j = 0; j < conf.num; ++j) {
+ if (conf.type[j] == i)
+ tmp[currId++] = j;
+ }
+
+ // Loop over rigid body grid names associated with particle type
+ const std::vector<String>& gridNames = conf.partRigidBodyGrid[i];
+ for (int j = 0; j < gridNames.size(); ++j) {
+
+ // Loop over RB types
+ for (int rb = 0; rb < conf.numRigidTypes; ++rb) {
+ RigidBodyType& t = conf.rigidBody[rb];
+ const std::vector<String>& keys = t.potentialGridKeys;
+
+ // Loop over potential grids
+ for(int k = 0; k < keys.size(); k++) {
+ // printf(" checking grid keys ");
+ if (gridNames[j] == keys[k]) {
+ memcpy( &(t.particles[k][t.numParticles[k]]), tmp, sizeof(int)*currId );
+ t.numParticles[k] += currId;
+ }
+ }
+ }
+ }
+ }
+
+ // Initialize device data for RB force pairs after std::vector is done growing
+
+ // for (int i = 0; i < forcePairs.size(); i++)
+ // forcePairs[i].initialize();
+
+}
+
+
+void RigidBodyController::updateForces(Vector3* pos_d, Vector3* force_d, int s) {
if (s <= 1)
gpuErrchk( cudaProfilerStart() );
@@ -181,6 +264,12 @@ void RigidBodyController::updateForces(int s) {
}
}
+ // Grid–particle forces
+ for (int i = 0; i < rigidBodyByType.size(); i++) {
+ callGridParticleForceKernel( pos_d, force_d, conf.rigidBody[i], rigidBodyByType[i], s );
+ }
+
+ // Grid–Grid forces
if (forcePairs.size() > 0) {
for (int i=0; i < forcePairs.size(); i++)
@@ -204,7 +293,7 @@ void RigidBodyController::updateForces(int s) {
}
}
void RigidBodyController::integrate(int step) {
- // tell RBs to integrate
+ // tell RBs to integrate
for (int i = 0; i < rigidBodyByType.size(); i++) {
for (int j = 0; j < rigidBodyByType[i].size(); j++) {
RigidBody& rb = rigidBodyByType[i][j];
@@ -353,6 +442,69 @@ void RigidBodyForcePair::callGridForceKernel(int pairId, int s) {
lastRbGridID = i;
}
}
+void RigidBodyController::callGridParticleForceKernel(Vector3* pos_d, Vector3* force_d,
+ const RigidBodyType& t, std::vector<RigidBody>& rbs, int s) {
+ // get the force/torque on a rigid body, and forces on particles
+
+ // RBTODO: consolidate CUDA stream management
+ for (int i = 0; i < t.numPotGrids; ++i) {
+ if (t.numParticles[i] == 0) continue;
+
+ for (int j = 0; j < rbs.size(); ++j) {
+ // const int nb = 500;
+ /*
+ r: postion of particle in real space
+ B: grid Basis
+ o: grid origin
+ R: rigid body orientation
+ c: rigid body center
+
+ B': R.B
+ c': R.o + c
+ */
+ // Matrix3 B1 = getBasis1(i);
+ Vector3 c = rbs[j].getOrientation()*t.potentialGrids[i].getOrigin() + rbs[j].getPosition();
+ Matrix3 B = (rbs[j].getOrientation()*t.potentialGrids[i].getBasis()).inverse();
+
+ // RBTODO: get energy
+ const int nb = (t.numParticles[i]/NUMTHREADS)+1;
+
+ // RBTODO: IMPORTANT: Improve this
+ Vector3 forces[nb];
+ Vector3 torques[nb];
+ for (int k=0; k < nb; ++k) {
+ forces[k] = Vector3(0.0f);
+ torques[k] = Vector3(0.0f);
+ }
+ Vector3* forces_d;
+ Vector3* torques_d;
+ gpuErrchk(cudaMalloc(&forces_d, sizeof(Vector3)*nb));
+ gpuErrchk(cudaMalloc(&torques_d, sizeof(Vector3)*nb));
+ gpuErrchk(cudaMemcpy(forces_d, forces, sizeof(Vector3)*nb, cudaMemcpyHostToDevice));
+ gpuErrchk(cudaMemcpy(torques_d, torques, sizeof(Vector3)*nb, cudaMemcpyHostToDevice));
+
+ computePartGridForce<<< nb, NUMTHREADS, NUMTHREADS*2*sizeof(Vector3) >>>(
+ pos_d, force_d, t.numParticles[i], t.particles[i],
+ t.rawPotentialGrids_d[i],
+ B, c, forces_d, torques_d);
+
+ gpuErrchk(cudaMemcpy(forces, forces_d, sizeof(Vector3)*nb, cudaMemcpyDeviceToHost));
+ gpuErrchk(cudaMemcpy(torques, torques_d, sizeof(Vector3)*nb, cudaMemcpyDeviceToHost));
+
+ Vector3 f = Vector3(0.0f);
+ Vector3 t = Vector3(0.0f);
+ for (int k = 0; k < nb; ++k) {
+ f = f + forces[k];
+ t = t + torques[j];
+ }
+
+ t = -t - (rbs[j].getPosition()-c).cross( -f );
+ rbs[j].addForce( -f );
+ rbs[j].addTorque( t );
+ }
+ }
+}
+
void RigidBodyForcePair::retrieveForcesForGrid(const int i) {
// i: grid ID (less than numGrids)
const cudaStream_t &s = stream[streamID[i]];
diff --git a/RigidBodyController.h b/RigidBodyController.h
index f04251f606ea6350abeffa54f2c10a9c2d4898cc..8aed324cd00d4ea7605a9264a00b625cfc44179e 100644
--- a/RigidBodyController.h
+++ b/RigidBodyController.h
@@ -94,11 +94,13 @@ public:
RigidBodyController(const Configuration& c, const char* outArg);
void integrate(int step);
- void updateForces(int s);
+ void updateForces(Vector3* pos_d, Vector3* force_d, int s);
+ void callGridParticleForceKernel(Vector3* pos_d, Vector3* force_d, const RigidBodyType& t, std::vector<RigidBody>& rbs, int s);
private:
void copyGridsToDevice();
void initializeForcePairs();
+ void initializeParticleLists();
void print(int step);
void printLegend(std::ofstream &file);
diff --git a/RigidBodyType.cu b/RigidBodyType.cu
index b8ea109855c6dc2cdfc80ddf98b4287acbcb232b..c3f600d9d0be96174900453ee3b465e1a29f9ed7 100644
--- a/RigidBodyType.cu
+++ b/RigidBodyType.cu
@@ -18,6 +18,14 @@ void RigidBodyType::clear() {
potentialGridKeys.clear();
densityGridKeys.clear();
pmfKeys.clear();
+
+ if (numParticles > 0) {
+ for (int i=0; i < numParticles[i]; ++i) {
+ delete[] particles[i];
+ }
+ delete[] numParticles;
+ delete[] particles;
+ }
if (numPotGrids > 0) delete[] rawPotentialGrids;
if (numDenGrids > 0) delete[] rawDensityGrids;
@@ -26,44 +34,6 @@ void RigidBodyType::clear() {
}
-// void RigidBodyType::copy(const RigidBodyType& src) {
-// this = new RigidBodyType(src.name);
-// num = src.num;
-// // if (src.pmf != NULL) pmf = new BaseGrid(*src.pmf);
-// if (src.reservoir != NULL) reservoir = new Reservoir(*src.reservoir);
-
-// numPotentialGrid = src.numPotentialGrid;
-// // TODO: fix this: BaseGrid*[]
-// for (int i=0; i < numPotentialGrid; i++) {
-
-// }
-
-// numDensityGrid = src.numDensityGrid;
-// }
-
-// RigidBodyType& RigidBodyType::operator=(const RigidBodyType& src) {
-// clear();
-// copy(src);
-// return *this;
-// }
-
-
-// KeyGrid RigidBodyType::createKeyGrid(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);
-// KeyGrid g;
-// g.key = token[0];
-// g.grid = * new BaseGrid(token[1]);
-// return g;
-// }
-
-
// void RigidBodyType::setDampingCoeffs(float timestep, float tmp_mass, Vector3 tmp_inertia, float tmp_transDamping, float tmp_rotDamping) {
void RigidBodyType::setDampingCoeffs(float timestep) { /* MUST ONLY BE CALLED ONCE!!! */
/*–––––––––––––––––––––––––––––––––––––––––––––––––––––––.
diff --git a/RigidBodyType.h b/RigidBodyType.h
index 44a813ded4ee4f6112a91f898b7ef88658a4eef7..4d9849dea3fe999636d92417581c55348a1cb724 100644
--- a/RigidBodyType.h
+++ b/RigidBodyType.h
@@ -40,7 +40,7 @@ RigidBodyType(const String& name = "") :
name(name), num(0),
reservoir(NULL), mass(1.0f), inertia(), transDamping(),
rotDamping(), initPos(), initRot(Matrix3(1.0f)),
- numPotGrids(0), numDenGrids(0), numPmfs(0) { }
+ numPotGrids(0), numDenGrids(0), numPmfs(0), numParticles(0) { }
/* RigidBodyType(const RigidBodyType& src) { copy(src); } */
~RigidBodyType() { clear(); }
@@ -56,7 +56,8 @@ RigidBodyType(const String& name = "") :
void updateRaw();
void setDampingCoeffs(float timestep);
-
+
+ // TODO: privatize
public:
@@ -91,9 +92,12 @@ public:
std::vector<float> densityGridScale;
std::vector<float> pmfScale;
+ int* numParticles; /* particles affected by potential grids */
+ int** particles;
// RBTODO: clear std::vectors after initialization, (but keep offsets)
// duplicates of std::vector grids for device
+public:
int numPotGrids;
int numDenGrids;
int numPmfs;
@@ -105,6 +109,7 @@ public:
Vector3* rawPotentialOrigins;
Vector3* rawDensityOrigins;
+
// device pointers
RigidBodyGrid** rawPotentialGrids_d;
RigidBodyGrid** rawDensityGrids_d;