diff --git a/src/BaseGrid.h b/src/BaseGrid.h
index a2c4d588b577df9fa9df796bbd559b82b464f9e1..55669bd770225d88c910ce94f3f6c9e7e15dc348 100644
--- a/src/BaseGrid.h
+++ b/src/BaseGrid.h
@@ -21,6 +21,18 @@
#include <ctime>
// #include <cuda.h>
+enum BoundaryCondition { dirichlet, neumann, periodic };
+enum InterpolationOrder { linear = 1, cubic = 3 };
+
+#define INTERPOLATE_FORCE(result, function, boundary_condition, args) \
+switch (boundary_condition) { \
+case dirichlet: \
+ result = function<dirichlet>(args); break; \
+case neumann: \
+ result = function<neumann>(args); break; \
+case periodic: \
+ result = function<periodic>(args); break; \
+}
// using namespace std;
@@ -500,6 +512,7 @@ public:
}
+ template <BoundaryCondition bc>
DEVICE inline ForceEnergy interpolateForceDLinearly(const Vector3& pos) const {
const Vector3 l = basisInv.transform(pos - origin);
@@ -512,6 +525,13 @@ public:
const float wy = l.y - homeY;
const float wz = l.z - homeZ;
+ if (bc == neumann) {
+ if (homeX < -1 || homeX >= nx+1 ||
+ homeY < -1 || homeY >= ny+1 ||
+ homeZ < -1 || homeZ >= nz+1)
+ return ForceEnergy();
+ }
+
float v[2][2][2];
for (int iz = 0; iz < 2; iz++) {
int jz = (iz + homeZ);
@@ -519,9 +539,31 @@ public:
int jy = (iy + homeY);
for (int ix = 0; ix < 2; ix++) {
int jx = (ix + homeX);
- int ind = jz + jy*nz + jx*nz*ny;
- v[ix][iy][iz] = jz < 0 || jz >= nz || jy < 0 || jy >= ny || jx < 0 || jx >= nx ?
+ int ind;
+ switch (bc) {
+ case dirichlet:
+ ind = jz + jy*nz + jx*nz*ny;
+ v[ix][iy][iz] =
+ jz < 0 || jz >= nz ||
+ jy < 0 || jy >= ny ||
+ jx < 0 || jx >= nx ?
0 : val[ind];
+ break;
+ case neumann:
+ ind =
+ (jz < 0 ? 0 : jz >= nz ? nz-1 : jz) +
+ (jy < 0 ? 0 : jy >= ny ? ny-1 : jy)*nz +
+ (jx < 0 ? 0 : jx >= nx ? nx-1 : jx)*nz*ny;
+ v[ix][iy][iz] = val[ind];
+ break;
+ case periodic:
+ ind =
+ (jz < 0 ? nz-1 : jz >= nz ? 0 : jz) +
+ (jy < 0 ? ny-1 : jy >= ny ? 0 : jy)*nz +
+ (jx < 0 ? nx-1 : jx >= nx ? 0 : jx)*nz*ny;
+ v[ix][iy][iz] = val[ind];
+ break;
+ }
}
}
}
@@ -865,71 +907,6 @@ public:
return val[j];*/
}
- //#define cubic
- DEVICE inline ForceEnergy interpolateForceDLinearlyPeriodic(const Vector3& pos) const {
- //#ifdef cubic
- //return interpolateForceD(pos);
- //#elif defined(cubic_namd)
- //return interpolateForceDnamd(pos);
- //#else
- return interpolateForceDLinearly(pos);
- //#endif
- #if 0
- const Vector3 l = basisInv.transform(pos - origin);
-
- // Find the home node.
- const int homeX = int(floor(l.x));
- const int homeY = int(floor(l.y));
- const int homeZ = int(floor(l.z));
-
- const float wx = l.x - homeX;
- const float wy = l.y - homeY;
- const float wz = l.z - homeZ;
-
- float v[2][2][2];
- for (int iz = 0; iz < 2; iz++) {
- int jz = (iz + homeZ);
- jz = (jz < 0) ? nz-1 : jz;
- jz = (jz >= nz) ? 0 : jz;
- for (int iy = 0; iy < 2; iy++) {
- int jy = (iy + homeY);
- jy = (jy < 0) ? ny-1 : jy;
- jy = (jy >= ny) ? 0 : jy;
- for (int ix = 0; ix < 2; ix++) {
- int jx = (ix + homeX);
- jx = (jx < 0) ? nx-1 : jx;
- jx = (jx >= nx) ? 0 : jx;
- int ind = jz + jy*nz + jx*nz*ny;
- v[ix][iy][iz] = val[ind];
- // printf("%d %d %d: %d %f\n",ix,iy,iz,ind,val[ind]); looks OK
- }
- }
- }
-
- float g3[3][2];
- for (int iz = 0; iz < 2; iz++) {
- float g2[2][2];
- for (int iy = 0; iy < 2; iy++) {
- g2[0][iy] = (v[1][iy][iz] - v[0][iy][iz]); /* f.x */
- g2[1][iy] = wx * (v[1][iy][iz] - v[0][iy][iz]) + v[0][iy][iz]; /* f.y & f.z */
- }
- // Mix along y.
- g3[0][iz] = wy * (g2[0][1] - g2[0][0]) + g2[0][0];
- g3[1][iz] = (g2[1][1] - g2[1][0]);
- g3[2][iz] = wy * (g2[1][1] - g2[1][0]) + g2[1][0];
- }
- // Mix along z.
- Vector3 f;
- f.x = -(wz * (g3[0][1] - g3[0][0]) + g3[0][0]);
- f.y = -(wz * (g3[1][1] - g3[1][0]) + g3[1][0]);
- f.z = - (g3[2][1] - g3[2][0]);
-
- f = basisInv.transpose().transform(f);
- float e = wz * (g3[2][1] - g3[2][0]) + g3[2][0];
- return ForceEnergy(f,e);
- #endif
- }
-
inline virtual Vector3 interpolateForce(Vector3 pos) const {
Vector3 f;
Vector3 l = basisInv.transform(pos - origin);
diff --git a/src/BrownianParticleType.cpp b/src/BrownianParticleType.cpp
index 50abb69740c2b6fdc4bad1a9083c7ae0425a943d..e514500ec3b6b87058cf58389c2a99ab7ddd506d 100644
--- a/src/BrownianParticleType.cpp
+++ b/src/BrownianParticleType.cpp
@@ -12,6 +12,7 @@ void BrownianParticleType::clear() {
if (reservoir != NULL) delete reservoir;
if (meanPmf != NULL) delete [] meanPmf;
pmf = NULL, diffusionGrid = NULL;
+ pmf_boundary_conditions = NULL;
forceXGrid = NULL, forceYGrid = NULL, forceZGrid = NULL;
reservoir = NULL, meanPmf = NULL;
}
@@ -25,6 +26,7 @@ void BrownianParticleType::copy(const BrownianParticleType& src) {
radius = src.radius;
eps = src.eps;
pmf = src.pmf;
+ pmf_boundary_conditions = src.pmf_boundary_conditions;
meanPmf = src.meanPmf;
numPartGridFiles = src.numPartGridFiles;
//Han-Yi Chou
diff --git a/src/BrownianParticleType.h b/src/BrownianParticleType.h
index 68f99a9bb03d7cec972931bbebcb9347fef3f852..524a5b66be452771f0aa8aec65d76dfddb60b7dd 100644
--- a/src/BrownianParticleType.h
+++ b/src/BrownianParticleType.h
@@ -30,7 +30,8 @@ class BrownianParticleType {
BrownianParticleType(const String& name = "") :
name(name), num(0),
diffusion(0.0f), radius(1.0f), charge(0.0f), eps(0.0f), meanPmf(NULL),
- numPartGridFiles(-1), reservoir(NULL), pmf(NULL), diffusionGrid(NULL),
+ numPartGridFiles(-1), reservoir(NULL), pmf(NULL), pmf_boundary_conditions(NULL),
+ diffusionGrid(NULL),
forceXGrid(NULL), forceYGrid(NULL), forceZGrid(NULL){ }
BrownianParticleType(const BrownianParticleType& src) { copy(src); }
@@ -62,6 +63,7 @@ public:
Reservoir* reservoir;
BaseGrid* pmf;
+ BoundaryCondition* pmf_boundary_conditions;
BaseGrid* diffusionGrid;
BaseGrid* forceXGrid;
BaseGrid* forceYGrid;
diff --git a/src/Configuration.cpp b/src/Configuration.cpp
index 9da0ba8cc449d9b5539a467ffd64a9599e4f4737..6e9f50194be10071d3f42263b891b374d7966433 100644
--- a/src/Configuration.cpp
+++ b/src/Configuration.cpp
@@ -599,18 +599,28 @@ void Configuration::copyToCUDA() {
for (int i = 0; i < numParts; i++)
{
- BaseGrid *pmf = NULL, *diffusionGrid = NULL;
+ BaseGrid *pmf = NULL, *diffusionGrid = NULL;
BrownianParticleType *b = new BrownianParticleType(part[i]);
// Copy PMF
if (part[i].pmf != NULL)
{
- float *tmp;
- gpuErrchk(cudaMalloc(&pmf, sizeof(BaseGrid)*part[i].numPartGridFiles));
- gpuErrchk(cudaMalloc(&tmp, sizeof(float)*part[i].numPartGridFiles));
- gpuErrchk(cudaMemcpy(tmp, part[i].meanPmf, sizeof(float)*part[i].numPartGridFiles,
- cudaMemcpyHostToDevice));
- b->meanPmf = tmp;
+ {
+ float *tmp;
+ gpuErrchk(cudaMalloc(&tmp, sizeof(float)*part[i].numPartGridFiles));
+ gpuErrchk(cudaMemcpy(tmp, part[i].meanPmf, sizeof(float)*part[i].numPartGridFiles,
+ cudaMemcpyHostToDevice));
+ b->meanPmf = tmp;
+ }
+
+ {
+ BoundaryCondition *tmp;
+ size_t s = sizeof(BoundaryCondition)*part[i].numPartGridFiles;
+ gpuErrchk(cudaMalloc(&tmp, s));
+ gpuErrchk(cudaMemcpy(tmp, part[i].pmf_boundary_conditions, s, cudaMemcpyHostToDevice));
+ b->pmf_boundary_conditions = tmp;
+ }
+ gpuErrchk(cudaMalloc(&pmf, sizeof(BaseGrid)*part[i].numPartGridFiles));
for(int j = 0; j < part[i].numPartGridFiles; ++j)
{
float *val = NULL;
@@ -626,6 +636,7 @@ void Configuration::copyToCUDA() {
}
b->pmf = pmf;
+
// Copy the diffusion grid
if (part[i].diffusionGrid != NULL) {
float *val = NULL;
@@ -1049,6 +1060,27 @@ int Configuration::readParameters(const char * config_file) {
//partGridFileScale[currPart] = (float) strtod(value.val(), NULL);
stringToArray<float>(&value, part[currPart].numPartGridFiles,
&partGridFileScale[currPart]);
+ } else if (param == String("gridFileBoundaryConditions")) {
+ register size_t num = value.tokenCount();
+ String *tokens = new String[num];
+ BoundaryCondition *data = new BoundaryCondition[num];
+ value.tokenize(tokens);
+ for(size_t i = 0; i < num; ++i) {
+ tokens[i].lower();
+ if (tokens[i] == "dirichlet")
+ data[i] = dirichlet;
+ else if (tokens[i] == "neumann")
+ data[i] = neumann;
+ else if (tokens[i] == "periodic")
+ data[i] = periodic;
+ else {
+ fprintf(stderr,"WARNING: Unrecognized gridFile boundary condition \"%s\". Using Dirichlet.\n", tokens[i].val() );
+ data[i] = dirichlet;
+ }
+ }
+ delete[] tokens;
+ part[currPart].pmf_boundary_conditions = data;
+
} else if (param == String("rigidBodyPotential")) {
partRigidBodyGrid[currPart].push_back(value);
}
@@ -1119,12 +1151,22 @@ int Configuration::readParameters(const char * config_file) {
//partGridFile[currPart] = value;
stringToArray<String>(&value, part[currPart].numPartGridFiles,
&partGridFile[currPart]);
- partGridFileScale[currPart] = new float[part[currPart].numPartGridFiles];
- for(int i = 0; i < part[currPart].numPartGridFiles; ++i) {
+ const int& num = part[currPart].numPartGridFiles;
+ partGridFileScale[currPart] = new float[num];
+ for(int i = 0; i < num; ++i) {
printf("%s ", partGridFile[currPart]->val());
partGridFileScale[currPart][i] = 1.0f;
}
+ // Set default boundary conditions for grids
+ BoundaryCondition *bc = part[currPart].pmf_boundary_conditions;
+ if (bc == NULL) {
+ bc = new BoundaryCondition[num];
+ for(int i = 0; i < num; ++i) {
+ bc[i] = dirichlet;
+ }
+ part[currPart].pmf_boundary_conditions = bc;
+ }
}
else if (currPartClass == partClassRB)
rigidBody[currRB].addPMF(value);
diff --git a/src/GrandBrownTown.cuh b/src/GrandBrownTown.cuh
index 2dc5d5c146ed5f282083ca9063ef0a4f858840ed..ba31f468596616e431b35c7bc82998ef69659c92 100644
--- a/src/GrandBrownTown.cuh
+++ b/src/GrandBrownTown.cuh
@@ -118,9 +118,10 @@ updateKernelNoseHooverLangevin(Vector3* __restrict__ pos, Vector3* __restrict__
for(int i = 0; i < pt.numPartGridFiles; ++i)
{
ForceEnergy tmp(0.f,0.f);
- if(!scheme)
- tmp = pt.pmf[i].interpolateForceDLinearlyPeriodic(r0);
- else
+ if(!scheme) {
+ const BoundaryCondition& bc = pt.pmf_boundary_conditions[i];
+ INTERPOLATE_FORCE(tmp, pt.pmf[i].interpolateForceDLinearly, bc, r0)
+ } else
tmp = pt.pmf[i].interpolateForceD(r0);
fe.f += tmp.f;
}
@@ -170,7 +171,7 @@ updateKernelNoseHooverLangevin(Vector3* __restrict__ pos, Vector3* __restrict__
p2 = sys->wrapDiff( p2 ) + gridCenter;
ForceEnergy diff;
if(!scheme)
- diff = pt.diffusionGrid->interpolateForceDLinearlyPeriodic(p2);
+ diff = pt.diffusionGrid->interpolateForceDLinearly<periodic>(p2);
else
diff = pt.diffusionGrid->interpolateForceD(p2);
gamma = Vector3(kTlocal / (mass * diff.e));
@@ -233,9 +234,10 @@ __global__ void updateKernelBAOAB(Vector3* pos, Vector3* momentum, Vector3* __re
for(int i = 0; i < pt.numPartGridFiles; ++i)
{
ForceEnergy tmp(0.f, 0.f);
- if(!scheme)
- tmp = pt.pmf[i].interpolateForceDLinearlyPeriodic(r0);
- else
+ if(!scheme) {
+ BoundaryCondition bc = pt.pmf_boundary_conditions[i];
+ INTERPOLATE_FORCE(tmp, pt.pmf[i].interpolateForceDLinearly, bc, r0)
+ } else
tmp = pt.pmf[i].interpolateForceD(r0);
fe.f += tmp.f;
}
@@ -282,7 +284,7 @@ __global__ void updateKernelBAOAB(Vector3* pos, Vector3* momentum, Vector3* __re
p2 = sys->wrapDiff( p2 ) + gridCenter;
ForceEnergy diff;
if(!scheme)
- diff = pt.diffusionGrid->interpolateForceDLinearlyPeriodic(p2);
+ diff = pt.diffusionGrid->interpolateForceDLinearly<periodic>(p2);
else
diff = pt.diffusionGrid->interpolateForceD(p2);
gamma = Vector3(kTlocal / (mass * diff.e));
@@ -332,9 +334,10 @@ __global__ void LastUpdateKernelBAOAB(Vector3* pos,Vector3* momentum, Vector3* _
for(int i = 0; i < pt.numPartGridFiles; ++i)
{
ForceEnergy tmp(0.f, 0.f);
- if(!scheme)
- tmp = pt.pmf[i].interpolateForceDLinearlyPeriodic(r0);
- else
+ if(!scheme) {
+ BoundaryCondition bc = pt.pmf_boundary_conditions[i];
+ INTERPOLATE_FORCE(tmp, pt.pmf[i].interpolateForceDLinearly, bc, r0)
+ } else
tmp = pt.pmf[i].interpolateForceD(r0);
fe += tmp;
@@ -463,9 +466,10 @@ void updateKernel(Vector3* pos, Vector3* __restrict__ forceInternal, int type[],
for(int i = 0; i < pt.numPartGridFiles; ++i)
{
ForceEnergy tmp(0.f,0.f);
- if(!scheme)
- tmp = pt.pmf[i].interpolateForceDLinearlyPeriodic(p);
- else
+ if(!scheme) {
+ BoundaryCondition bc = pt.pmf_boundary_conditions[i];
+ INTERPOLATE_FORCE(tmp, pt.pmf[i].interpolateForceDLinearly, bc, p)
+ } else
tmp = pt.pmf[i].interpolateForceD(p);
fe += tmp;
}
@@ -525,7 +529,7 @@ void updateKernel(Vector3* pos, Vector3* __restrict__ forceInternal, int type[],
/* printf("atom %d: ps2: %f %f %f\n", idx, p2.x, p2.y, p2.z); */
ForceEnergy diff;
if(!scheme)
- diff = pt.diffusionGrid->interpolateForceDLinearlyPeriodic(p2);
+ diff = pt.diffusionGrid->interpolateForceDLinearly<periodic>(p2);
else
diff = pt.diffusionGrid->interpolateForceD(p2);
diffusion = diff.e;