From ad21c2ddb05adc5cc9c8313e3be43c551f8bf615 Mon Sep 17 00:00:00 2001
From: Emmanuel Guzman <eguzman6@illinois.edu>
Date: Thu, 30 Jun 2016 14:24:24 -0500
Subject: [PATCH] computeTabulateBonds works, bug computing forces fixed,
bondList is created, but not implemented yet.
---
ComputeForce.cu | 4 ++--
ComputeForce.cuh | 32 ++++++++++++++++++++-----------
Configuration.h | 3 ++-
GrandBrownTown.cu | 45 ++++++++++++++++++++++++++++++++++++++++++--
GrandBrownTown.h | 4 +++-
TabulatedPotential.h | 6 +++---
makefile | 4 ++--
7 files changed, 76 insertions(+), 22 deletions(-)
diff --git a/ComputeForce.cu b/ComputeForce.cu
index 7dce059..81e6726 100644
--- a/ComputeForce.cu
+++ b/ComputeForce.cu
@@ -591,10 +591,10 @@ float ComputeForce::computeTabulated(Vector3* force, Vector3* pos, int* type,
tableAngle_d, numAngles, num, sys_d, energies_d, get_energy);
/***************************************** call computeTabulatedBonds *****************************************/
+ //if(bondMap != NULL && tableBond_d != NULL)
if(bondMap != NULL && tableBond_d != NULL)
{
- computeTabulatedBonds <<<numBlocks, numThreads>>> ( force, pos, num, numParts, sys_d, bonds, bondMap, numBonds, numReplicas, energies_d, get_energy, tableBond);
-
+ computeTabulatedBonds <<<numBlocks, numThreads>>> ( force, pos, num, numParts, sys_d, bonds, bondMap, numBonds, numReplicas, energies_d, get_energy, tableBond_d);
}
/***************************************** end *****************************************/
computeDihedrals<<<numBlocks, numThreads>>>(force, pos, dihedrals,
diff --git a/ComputeForce.cuh b/ComputeForce.cuh
index 4362878..debebd7 100644
--- a/ComputeForce.cuh
+++ b/ComputeForce.cuh
@@ -449,7 +449,7 @@ __global__ void computeTabulatedKernel(
// Particle's type and position
const int ind = g_pairTabPotType[i];
- /* printf("tid,bid,pos[ai(%d)]: %d %d %f %f %f\n", ai, threadIdx.x, blockIdx.x, pos[ai].x, pos[ai].y, pos[ai].z); */
+ /* printf("tid,bid,pos[ai(%d)]: %d %d %f %f %f\n", ai, threadIdx.x, blockIdx.x, pos[ai].x, pos[ai].y, pos[ai].z); //*/
Vector3 dr = pos[aj] - pos[ai];
dr = sys->d_wrapDiff(dr);
@@ -458,8 +458,8 @@ __global__ void computeTabulatedKernel(
float d2 = dr.length2();
if (tablePot[ind] != NULL && d2 <= cutoff2) {
Vector3 f = tablePot[ind]->computef(dr,d2);
- atomicAdd( &force[ai], -f );
- atomicAdd( &force[aj], f );
+ atomicAdd( &force[ai], f );
+ atomicAdd( &force[aj], -f );
}
}
}
@@ -646,7 +646,7 @@ void computeAngles(Vector3 force[], Vector3 pos[],
/***************************************** computeBonds *****************************************/
__global__ void computeTabulatedBonds( Vector3* force, Vector3* pos, int num,
- int numParts, BaseGrid* sys, Bond* bonds,
+ int numParts, BaseGrid* sys, Bond bonds[],
int2* bondMap, int numBonds, int numReplicas,
float* g_energies, bool get_energy, TabulatedPotential** tableBond)
{
@@ -659,7 +659,7 @@ __global__ void computeTabulatedBonds( Vector3* force, Vector3* pos, int num,
// Initialize interaction energy (per particle)
float energy_local = 0.0f;
- const int repID = i / num;
+ const int repID = i / num; //TODO: ask why repID is even here, it's unnecessary in this kernel
// BONDS
/* Each particle may have a varying number of bonds
@@ -698,18 +698,28 @@ __global__ void computeTabulatedBonds( Vector3* force, Vector3* pos, int num,
// tabulated potential
// If the user has specified the ADD option, then add the bond
// force to the tabulated potential value
- EnergyForce ft = tableBond[ bonds[currBond].tabFileIndex ]->compute(dr);
-
+
+ /* start testing */
+ //cout << "tableBond points to " << tableBond << ".\n"
+ //<< "tableBond has a size of " << sizeof(tableBond)<< ".\n";
+ //Bond test_bond = bonds[currBond];
+ //int test_num = test_bond.tabFileIndex;
+ //TabulatedPotential *test_pot = tableBond[test_num];
+ //EnergyForce test_energy = test_pot -> compute(dr);
+ /* end testing */
+
+ force_local += tableBond[ bonds[currBond].tabFileIndex ]->computef(dr,dr.length2());
+
+ /*EnergyForce ft = tableBond[ bonds[currBond].tabFileIndex ]->compute(dr);
force_local += ft.f;
if (get_energy && j > i)
- energy_local += ft.e;
+ energy_local += ft.e;*/
}
-
- force[i] = force_local;
+ atomicAdd( &force[i], force_local );
if (get_energy)
- g_energies[i] = energy_local;
+ atomicAdd( &g_energies[i], energy_local);
}
}
/***************************************** end *****************************************/
diff --git a/Configuration.h b/Configuration.h
index b9143b7..4829e40 100644
--- a/Configuration.h
+++ b/Configuration.h
@@ -109,6 +109,7 @@ public:
float timeLast; // used with posLast
float minimumSep; // minimum separation allowed with placing new particles
+
// RigidBody variables
/* int numRB; */
/* std::vector< std::vector<RigidBody> > rbs; */
@@ -178,7 +179,7 @@ public:
String* bondTableFile;
int numTabBondFiles;
int2* bondMap;
-
+
Exclude* excludes;
int2* excludeMap;
String excludeRule;
diff --git a/GrandBrownTown.cu b/GrandBrownTown.cu
index aa062e6..e5700e1 100644
--- a/GrandBrownTown.cu
+++ b/GrandBrownTown.cu
@@ -124,6 +124,7 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
numTabBondFiles = c.numTabBondFiles;
bondMap = c.bondMap;
+ // TODO: bondList = c.bondList;
excludes = c.excludes;
excludeMap = c.excludeMap;
@@ -179,6 +180,7 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
fullLongRange, numBonds, numTabBondFiles, numExcludes, numAngles, numTabAngleFiles,
numDihedrals, numTabDihedralFiles, numReplicas);
+ // TODO: check for duplicate potentials
if (c.tabulatedPotential) {
printf("Loading the tabulated potentials...\n");
for (int p = 0; p < numParts*numParts; p++) {
@@ -216,6 +218,30 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
internal->addDihedralPotential(dihedralTableFile[p].val(), p, dihedrals, dihedrals_d);
}
+/***************************************** init. bondList *****************************************/
+ bondList = new int3[numBonds/2];
+ int j = 0;
+
+ for(int i = 0; i < numBonds; ++i)
+ {
+ if(bonds[i].ind1 < bonds[i].ind2)
+ {
+ bondList[j] = make_int3( bonds[i].ind1, bonds[i].ind2, bonds[i].tabFileIndex );
+ cout << "Displaying: bondList["<< j <<"].x = " << bondList[j].x << ".\n"
+ << "Displaying: bondList["<< j <<"].y = " << bondList[j].y << ".\n"
+ << "Displaying: bondList["<< j <<"].z = " << bondList[j].z << ".\n";
+ j++;
+ }
+ }
+
+ //TODO: memCopy bondList here to device, pass in list to decompose and change the function signitures. Afterwards edit the function to make use of the list.
+ createBondList();
+
+ // TODO: copy data to device (not here)
+ // TODO: use bondList in computeTabulatedBonds kernel
+
+/***************************************** end *****************************************/
+
forceInternal = new Vector3[num * numReplicas];
if (fullLongRange != 0)
@@ -237,7 +263,7 @@ GrandBrownTown::~GrandBrownTown() {
delete[] pos;
delete[] type;
delete[] serial;
-
+ delete[] bondList;
delete randoGen;
// Auxillary objects
@@ -248,6 +274,7 @@ GrandBrownTown::~GrandBrownTown() {
gpuErrchk(cudaFree(pos_d));
gpuErrchk(cudaFree(forceInternal_d));
gpuErrchk(cudaFree(randoGen_d));
+ gpuErrchk( cudaFree(bondList_d) );
}
// Run the Brownian Dynamics steps.
@@ -303,7 +330,6 @@ void GrandBrownTown::run() {
rt_timer_start(timer0);
rt_timer_start(timerS);
-
// We haven't done any steps yet.
// Do decomposition if we have to
if (fullLongRange == 0)
@@ -959,3 +985,18 @@ void GrandBrownTown::copyToCUDA() {
cudaMemcpyHostToDevice));
gpuErrchk(cudaDeviceSynchronize());
}
+
+void GrandBrownTown::createBondList()
+{
+ size_t size = (numBonds/2) * sizeof(int3);
+ gpuErrchk( cudaMalloc( &bondList_d, sizeof(size) ) );
+ gpuErrchk( cudaMemcpyAsync( bondList_d, bondList, sizeof(size), cudaMemcpyHostToDevice) );
+
+ for(int i = 0 ; i < numBonds/2 ; i++)
+ {
+ cout << "Displaying: bondList_d["<< i <<"].x = " << bondList[i].x << ".\n"
+ << "Displaying: bondList_d["<< i <<"].y = " << bondList[i].y << ".\n"
+ << "Displaying: bondList_d["<< i <<"].z = " << bondList[i].z << ".\n";
+
+ }
+}
diff --git a/GrandBrownTown.h b/GrandBrownTown.h
index 43ba4d7..dc8c3be 100644
--- a/GrandBrownTown.h
+++ b/GrandBrownTown.h
@@ -98,6 +98,7 @@ private:
void getDebugForce();
void copyToCUDA();
+ void createBondList();
public:
// Compute the current in nanoamperes.
@@ -153,7 +154,6 @@ private:
RigidBodyController RBC;
Vector3* rbPos; // rigid body positions
-
// CUDA device variables
Vector3 *pos_d, *forceInternal_d, *force_d;
int *type_d;
@@ -166,6 +166,7 @@ private:
int2* excludeMap_d;
Angle* angles_d;
Dihedral* dihedrals_d;
+ int3 *bondList_d;
// System parameters
String outputName;
@@ -232,6 +233,7 @@ private:
String* bondTableFile;
int numTabBondFiles;
int2* bondMap;
+ int3 *bondList;
Exclude* excludes;
int2* excludeMap;
diff --git a/TabulatedPotential.h b/TabulatedPotential.h
index 91a401a..1df4d69 100644
--- a/TabulatedPotential.h
+++ b/TabulatedPotential.h
@@ -77,7 +77,7 @@ public:
// Interpolate.
float energy = du*w+u0;
- Vector3 force = (-du/(d*dr))*r;
+ Vector3 force = (du/(d*dr))*r;
return EnergyForce(energy,force);
}
HOST DEVICE inline EnergyForce compute(Vector3 r, float d) {
@@ -95,7 +95,7 @@ public:
// Interpolate.
float energy = du*w+u0;
- Vector3 force = (-du/(d*dr))*r;
+ Vector3 force = (du/(d*dr))*r;
return EnergyForce(energy,force);
}
HOST DEVICE inline Vector3 computef(Vector3 r, float d) {
@@ -110,7 +110,7 @@ public:
if (home >= n) return Vector3(0.0f);
if (home+1 < n)
- return (-(v0[home+1]-v0[home])/(d*dr))*r;
+ return ((v0[home+1]-v0[home])/(d*dr))*r;
else
return Vector3(0.0f);
}
diff --git a/makefile b/makefile
index 8a9a05e..881f152 100644
--- a/makefile
+++ b/makefile
@@ -12,7 +12,7 @@ include ./findcudalib.mk
INCLUDE = $(CUDA_PATH)/include
-DEBUG = -g -O0
+# DEBUG = -g -G -O0
CC_FLAGS = -Wall -Wno-write-strings -I$(INCLUDE) $(DEBUG) -std=c++0x -pedantic
# NV_FLAGS = --maxrregcount 63 -Xptxas -v # -v,-abi=no
NV_FLAGS = -Xptxas -v # -v,-abi=no
@@ -50,7 +50,7 @@ CODE_20 := -arch=sm_20
# NV_FLAGS += $(CODE_10) $(CODE_12) $(CODE_20) $(CODE_30) $(CODE_35)
# NV_FLAGS += -arch=sm_35
# SM=52
-SM ?= 35
+SM ?= 30
NV_FLAGS += -gencode arch=compute_$(SM),code=compute_$(SM)
--
GitLab