Adding leanmd simplified as an example
authorNikhil Jain <nikhil@illinois.edu>
Mon, 24 Sep 2012 19:13:40 +0000 (14:13 -0500)
committerNikhil Jain <nikhil@illinois.edu>
Mon, 24 Sep 2012 19:13:40 +0000 (14:13 -0500)
examples/charm++/leanmd/Cell.cc [new file with mode: 0644]
examples/charm++/leanmd/Cell.h [new file with mode: 0644]
examples/charm++/leanmd/Compute.cc [new file with mode: 0644]
examples/charm++/leanmd/Compute.h [new file with mode: 0644]
examples/charm++/leanmd/Main.cc [new file with mode: 0644]
examples/charm++/leanmd/Main.h [new file with mode: 0644]
examples/charm++/leanmd/Makefile [new file with mode: 0644]
examples/charm++/leanmd/README [new file with mode: 0644]
examples/charm++/leanmd/defs.h [new file with mode: 0644]
examples/charm++/leanmd/leanmd.ci [new file with mode: 0644]
examples/charm++/leanmd/physics.h [new file with mode: 0644]

diff --git a/examples/charm++/leanmd/Cell.cc b/examples/charm++/leanmd/Cell.cc
new file mode 100644 (file)
index 0000000..d04389d
--- /dev/null
@@ -0,0 +1,280 @@
+
+#include "defs.h"
+#include "leanmd.decl.h"
+#include "Cell.h"
+
+extern /* readonly */ CProxy_Main mainProxy;
+extern /* readonly */ CProxy_Cell cellArray;
+extern /* readonly */ CProxy_Compute computeArray;
+extern /* readonly */ CkGroupID mCastGrpID;
+
+extern /* readonly */ int cellArrayDimX;
+extern /* readonly */ int cellArrayDimY;
+extern /* readonly */ int cellArrayDimZ;
+extern /* readonly */ int finalStepCount; 
+
+//default constructor
+Cell::Cell() {
+  __sdag_init();
+  int i;
+  inbrs = NUM_NEIGHBORS;
+  //load balancing to be called when AtSync is called
+  usesAtSync = CmiTrue;
+
+  int myid = thisIndex.x + thisIndex.y*cellArrayDimX + thisIndex.z*cellArrayDimX*cellArrayDimY;
+  myNumParts = PARTICLES_PER_CELL_START + (myid*(PARTICLES_PER_CELL_END-PARTICLES_PER_CELL_START))/(cellArrayDimX*cellArrayDimY*cellArrayDimZ);
+
+  // starting random generator
+  srand48(thisIndex.x+cellArrayDimX*(thisIndex.y+thisIndex.z*cellArrayDimY));
+
+  // Particle initialization
+  for(i=0; i < myNumParts; i++) {
+    particles.push_back(Particle());
+    particles[i].mass = HYDROGEN_MASS;
+
+    //give random values for position and velocity
+    particles[i].pos.x = drand48() * CELL_SIZE_X + thisIndex.x * CELL_SIZE_X;
+    particles[i].pos.y = drand48() * CELL_SIZE_Y + thisIndex.y * CELL_SIZE_Y;
+    particles[i].pos.z = drand48() * CELL_SIZE_Z + thisIndex.z * CELL_SIZE_Z;
+    particles[i].vel.x = (drand48() - 0.5) * .2 * MAX_VELOCITY;
+    particles[i].vel.y = (drand48() - 0.5) * .2 * MAX_VELOCITY;
+    particles[i].vel.z = (drand48() - 0.5) * .2 * MAX_VELOCITY;
+       particles[i].force.x = 0.0;
+       particles[i].force.y = 0.0;
+       particles[i].force.z = 0.0;
+  }
+
+  stepCount = 1;
+  updateCount = 0;
+  forceCount = 0;
+}
+
+//constructor for chare object migration
+Cell::Cell(CkMigrateMessage *msg): CBase_Cell(msg) {
+  __sdag_init();
+  usesAtSync = CmiTrue;
+  delete msg;
+}  
+
+Cell::~Cell() {}
+
+//function to create my computes
+void Cell::createComputes() {
+  int num;  
+
+  int x = thisIndex.x;
+  int y = thisIndex.y;
+  int z = thisIndex.z;
+  int px1, py1, pz1, dx, dy, dz, px2, py2, pz2;
+
+  // For Round Robin insertion
+  int numPes = CkNumPes();
+  int currPe = CkMyPe();
+
+  computesList = new int*[inbrs];
+  for (int i =0; i < inbrs; i++){
+    computesList[i] = new int[6];
+  }
+
+  /*  The computes X are inserted by a given cell:
+   *
+   *   ^  X  X  X
+   *   |  0  X  X
+   *   y  0  0  0
+   *      x ---->
+   */
+
+  for (num=0; num<inbrs; num++) {
+    dx = num / (NBRS_Y * NBRS_Z)            - NBRS_X/2;
+    dy = (num % (NBRS_Y * NBRS_Z)) / NBRS_Z   - NBRS_Y/2;
+    dz = num % NBRS_Z                       - NBRS_Z/2;
+
+    if (num >= inbrs/2){
+      px1 = x + 2;
+      py1 = y + 2;
+      pz1 = z + 2;
+      px2 = px1+dx;
+      py2 = py1+dy;
+      pz2 = pz1+dz;
+      computeArray(px1, py1, pz1, px2, py2, pz2).insert((++currPe)%numPes);
+      computesList[num][0] = px1; computesList[num][1] = py1; computesList[num][2] = pz1; 
+      computesList[num][3] = px2; computesList[num][4] = py2; computesList[num][5] = pz2;
+    }
+    else {
+      // these computes will be created by pairing cells
+      px1 = WRAP_X(x+dx)+2 ;
+      py1 = WRAP_Y(y+dy)+2;
+      pz1 = WRAP_Z(z+dz)+2;
+      px2 = px1 - dx;
+      py2 = py1 - dy;
+      pz2 = pz1 - dz;
+      computesList[num][0] = px1; computesList[num][1] = py1; computesList[num][2] = pz1; 
+      computesList[num][3] = px2; computesList[num][4] = py2; computesList[num][5] = pz2;
+    }
+  } // end of for loop
+  contribute(0,NULL,CkReduction::nop,CkCallback(CkReductionTarget(Main,computesCreated),mainProxy));
+}
+
+// Function to start interaction among particles in neighboring cells as well as its own particles
+void Cell::sendPositions() {
+  int len = particles.length();
+
+  //create the particle and control message to be sent to computes
+  ParticleDataMsg* msg = new (len) ParticleDataMsg;
+  msg->x = thisIndex.x;
+  msg->y = thisIndex.y;
+  msg->z = thisIndex.z;
+  msg->lengthAll = len;
+
+  for (int i = 0; i < len; i++)
+    msg->part[i] = particles[i].pos;
+
+  CkSetRefNum(msg,stepCount);
+
+  int px1, py1, pz1, px2, py2, pz2;
+
+  for(int num=0; num<inbrs; num++) {
+    px1 = computesList[num][0];
+    py1 = computesList[num][1];
+    pz1 = computesList[num][2];
+    px2 = computesList[num][3];
+    py2 = computesList[num][4];
+    pz2 = computesList[num][5];
+    if (num == inbrs-1)
+      computeArray(px1, py1, pz1, px2, py2, pz2).calculateForces(msg);
+    else {
+      ParticleDataMsg* newMsg = new (len) ParticleDataMsg;
+      newMsg->x = thisIndex.x;
+      newMsg->y = thisIndex.y;
+      newMsg->z = thisIndex.z;
+      newMsg->lengthAll = len;
+      memcpy(newMsg->part, msg->part, len*sizeof(vec3));
+      CkSetRefNum(newMsg,stepCount);
+      computeArray(px1, py1, pz1, px2, py2, pz2).calculateForces(newMsg);
+    }
+  }
+
+}
+
+//send the atoms that have moved beyond my cell to neighbors
+void Cell::migrateParticles(){
+  int i, x1, y1, z1;
+  CkVec<Particle> *outgoing = new CkVec<Particle>[inbrs];
+
+  for(i=0; i<particles.length(); i++) {
+    migrateToCell(particles[i], x1, y1, z1);
+    if(x1!=0 || y1!=0 || z1!=0) {
+      outgoing[(x1+1)*NBRS_Y*NBRS_Z + (y1+1)*NBRS_Z + (z1+1)].push_back(wrapAround(particles[i]));
+      particles.remove(i);
+    }
+  }
+  for(int num=0; num<inbrs; num++) {
+    x1 = num / (NBRS_Y * NBRS_Z)            - NBRS_X/2;
+    y1 = (num % (NBRS_Y * NBRS_Z)) / NBRS_Z - NBRS_Y/2;
+    z1 = num % NBRS_Z                       - NBRS_Z/2;
+    cellArray(WRAP_X(thisIndex.x+x1), WRAP_Y(thisIndex.y+y1), WRAP_Z(thisIndex.z+z1)).receiveParticles(outgoing[num]);
+  }
+  delete [] outgoing;
+}
+
+//check if the particle is to be moved
+void Cell::migrateToCell(Particle p, int &px, int &py, int &pz) {
+  // currently this is assuming that particles are
+  // migrating only to the immediate neighbors
+  int x = thisIndex.x * CELL_SIZE_X + CELL_ORIGIN_X;
+  int y = thisIndex.y * CELL_SIZE_Y + CELL_ORIGIN_Y;
+  int z = thisIndex.z * CELL_SIZE_Z + CELL_ORIGIN_Z;
+  px = py = pz = 0;
+
+  if (p.pos.x < x) px = -1;
+  else if (p.pos.x > x+CELL_SIZE_X) px = 1;
+
+  if (p.pos.y < y) py = -1;
+  else if (p.pos.y > y+CELL_SIZE_Y) py = 1;
+
+  if (p.pos.z < z) pz = -1;
+  else if (p.pos.z > z+CELL_SIZE_Z) pz = 1;
+}
+
+// Function to update properties (i.e. acceleration, velocity and position) in particles
+void Cell::updateProperties() {
+  int i;
+  double powTen, powFteen, realTimeDelta, invMassParticle;
+  powTen = pow(10.0, -10);
+  powFteen = pow(10.0, -15);
+  realTimeDelta = DEFAULT_DELTA * powFteen;
+  for(i = 0; i < particles.length(); i++) {
+    // applying kinetic equations
+    invMassParticle = 1 / particles[i].mass;
+    particles[i].acc = particles[i].force * invMassParticle;
+    particles[i].vel += particles[i].acc * realTimeDelta;
+
+    limitVelocity(particles[i]);
+
+    particles[i].pos += particles[i].vel * realTimeDelta;
+       particles[i].force.x = 0.0;
+       particles[i].force.y = 0.0;
+       particles[i].force.z = 0.0;
+  }
+  forceCount = 0;
+}
+
+void Cell::addForces(vec3 *forces){
+       // updating force information
+       for(int i = 0; i < particles.length(); i++){
+               particles[i].force += forces[i];
+       }
+}
+
+inline double velocityCheck(double inVelocity) {
+  if(fabs(inVelocity) > MAX_VELOCITY) {
+    if(inVelocity < 0.0 )
+      return -MAX_VELOCITY;
+    else
+      return MAX_VELOCITY;
+  } else {
+    return inVelocity;
+  }
+}
+
+void Cell::limitVelocity(Particle &p) {
+  p.vel.x = velocityCheck(p.vel.x);
+  p.vel.y = velocityCheck(p.vel.y);
+  p.vel.z = velocityCheck(p.vel.z);
+}
+
+Particle& Cell::wrapAround(Particle &p) {
+  if(p.pos.x < CELL_ORIGIN_X) p.pos.x += CELL_SIZE_X*cellArrayDimX;
+  if(p.pos.y < CELL_ORIGIN_Y) p.pos.y += CELL_SIZE_Y*cellArrayDimY;
+  if(p.pos.z < CELL_ORIGIN_Z) p.pos.z += CELL_SIZE_Z*cellArrayDimZ;
+  if(p.pos.x > CELL_ORIGIN_X + CELL_SIZE_X*cellArrayDimX) p.pos.x -= CELL_SIZE_X*cellArrayDimX;
+  if(p.pos.y > CELL_ORIGIN_Y + CELL_SIZE_Y*cellArrayDimY) p.pos.y -= CELL_SIZE_Y*cellArrayDimY;
+  if(p.pos.z > CELL_ORIGIN_Z + CELL_SIZE_Z*cellArrayDimZ) p.pos.z -= CELL_SIZE_Z*cellArrayDimZ;
+
+  return p;
+}
+
+//pack important data when I move
+void Cell::pup(PUP::er &p) {
+  CBase_Cell::pup(p);
+  __sdag_pup(p);
+  p | particles;
+  p | stepCount;
+  p | myNumParts;
+  p | updateCount;
+  p | inbrs;
+  p | forceCount;
+
+  if (p.isUnpacking()){
+    computesList = new int*[inbrs];
+    for (int i = 0; i < inbrs; i++){
+      computesList[i] = new int[6];
+    }
+  }
+
+  for (int i = 0; i < inbrs; i++){
+    PUParray(p, computesList[i], 6);
+  }
+
+}
+
diff --git a/examples/charm++/leanmd/Cell.h b/examples/charm++/leanmd/Cell.h
new file mode 100644 (file)
index 0000000..addf759
--- /dev/null
@@ -0,0 +1,58 @@
+#ifndef __CELL_H__
+#define __CELL_H__
+
+extern /* readonly */ CProxy_Main mainProxy;
+extern /* readonly */ CProxy_Cell cellArray;
+extern /* readonly */ CkGroupID mCastGrpID;
+extern /* readonly */ int finalStepCount;
+
+
+//data message to be sent to computes
+struct ParticleDataMsg : public CkMcastBaseMsg, public CMessage_ParticleDataMsg {
+  vec3* part; //list of atoms
+  int lengthAll;  //length of list
+  int x;    //x coordinate of cell sending this message
+  int y;    //y coordinate
+  int z;    //z coordinate
+
+  void pup(PUP::er &p){
+    CMessage_ParticleDataMsg::pup(p);
+    p | lengthAll;
+    p | x; p | y; p | z;
+    if (p.isUnpacking()){
+      part = new vec3[lengthAll];
+    }
+    PUParray(p, part, lengthAll);
+  } 
+};
+
+//chare used to represent a cell
+class Cell : public CBase_Cell {
+  private:
+    Cell_SDAG_CODE   //SDAG code
+    CkVec<Particle> particles;  //list of atoms
+    int **computesList;   //my compute locations
+    int stepCount;             // to count the number of steps, and decide when to stop
+    int myNumParts;   //number of atoms in my cell
+    int inbrs;        //number of interacting neighbors
+    int updateCount;
+       int forceCount;
+
+    void migrateToCell(Particle p, int &px, int &py, int &pz);
+    void updateProperties();   //updates properties after receiving forces from computes
+    void limitVelocity(Particle &p); //limit velcities to an upper limit
+    Particle& wrapAround(Particle &p); //particles going out of right enters from left
+    CProxySection_Compute mCastSecProxy; //handle to section proxy of computes
+
+  public:
+    Cell();
+    Cell(CkMigrateMessage *msg);
+    ~Cell();
+    void pup(PUP::er &p);
+    void createComputes();  //add my computes
+    void migrateParticles();
+    void sendPositions();
+       void addForces(vec3* forces);
+};
+
+#endif
diff --git a/examples/charm++/leanmd/Compute.cc b/examples/charm++/leanmd/Compute.cc
new file mode 100644 (file)
index 0000000..94286da
--- /dev/null
@@ -0,0 +1,45 @@
+#include "defs.h"
+#include "leanmd.decl.h"
+#include "Cell.h"
+#include "Compute.h"
+#include "physics.h"
+#include <algorithm>
+using std::swap;
+
+extern /* readonly */ CProxy_Main mainProxy;
+extern /* readonly */ CkGroupID mCastGrpID;
+
+extern /* readonly */ int cellArrayDimX;
+extern /* readonly */ int cellArrayDimY;
+extern /* readonly */ int cellArrayDimZ;
+extern /* readonly */ int finalStepCount; 
+
+//compute - Default constructor
+Compute::Compute() {
+  __sdag_init();
+  stepCount = 1;
+  usesAtSync = CmiTrue;
+}
+
+Compute::Compute(CkMigrateMessage *msg): CBase_Compute(msg)  { 
+  __sdag_init();
+  usesAtSync = CmiTrue;
+  delete msg;
+}
+
+//interaction within a cell
+void Compute::selfInteract(ParticleDataMsg *msg){
+  calcInternalForces(msg, stepCount);
+}
+
+//interaction between two cells
+void Compute::interact(ParticleDataMsg *msg){
+  calcPairForces(msg, bufferedMsg, stepCount);
+}
+
+//pack important information if I am moving
+void Compute::pup(PUP::er &p) {
+  CBase_Compute::pup(p);
+  __sdag_pup(p);
+  p | stepCount;
+}
diff --git a/examples/charm++/leanmd/Compute.h b/examples/charm++/leanmd/Compute.h
new file mode 100644 (file)
index 0000000..35bdd69
--- /dev/null
@@ -0,0 +1,22 @@
+#ifndef __COMPUTE_H__
+#define __COMPUTE_H__
+
+#include "defs.h"
+
+//class representing the interaction agents between a couple of cells
+class Compute : public CBase_Compute {
+  private:
+    Compute_SDAG_CODE
+    int stepCount;  //current step number
+    ParticleDataMsg *bufferedMsg; //copy of first message received for interaction
+
+  public:
+    Compute();
+    Compute(CkMigrateMessage *msg);
+    void pup(PUP::er &p);
+
+    void selfInteract(ParticleDataMsg *msg);
+    void interact(ParticleDataMsg *msg);
+};
+
+#endif
diff --git a/examples/charm++/leanmd/Main.cc b/examples/charm++/leanmd/Main.cc
new file mode 100644 (file)
index 0000000..5430547
--- /dev/null
@@ -0,0 +1,77 @@
+#include "time.h"
+
+#include "defs.h"
+#include "leanmd.decl.h"
+#include "Main.h"
+#include "Cell.h"
+#include "Compute.h"
+
+/* readonly */ CProxy_Main mainProxy;
+/* readonly */ CProxy_Cell cellArray;
+/* readonly */ CProxy_Compute computeArray;
+
+/* readonly */ int cellArrayDimX;
+/* readonly */ int cellArrayDimY;
+/* readonly */ int cellArrayDimZ;
+/* readonly */ int finalStepCount; 
+
+// Entry point of Charm++ application
+Main::Main(CkArgMsg* m) {
+   __sdag_init();
+  CkPrintf("\nLENNARD JONES MOLECULAR DYNAMICS START UP ...\n");
+
+  //set variable values to a default set
+  cellArrayDimX = CELLARRAY_DIM_X;
+  cellArrayDimY = CELLARRAY_DIM_Y;
+  cellArrayDimZ = CELLARRAY_DIM_Z;
+
+  mainProxy = thisProxy;
+
+  int numPes = CkNumPes();
+  int currPe = -1, pe;
+  int cur_arg = 1;
+
+  CkPrintf("\nInput Parameters...\n");
+
+  //read user parameters
+  //number of cells in each dimension
+  if (m->argc > cur_arg) {
+    cellArrayDimX=atoi(m->argv[cur_arg++]);
+    cellArrayDimY=atoi(m->argv[cur_arg++]);
+    cellArrayDimZ=atoi(m->argv[cur_arg++]);
+    CkPrintf("Cell Array Dimension X:%d Y:%d Z:%d of size %d %d %d\n",cellArrayDimX,cellArrayDimY,cellArrayDimZ,CELL_SIZE_X,CELL_SIZE_Y,CELL_SIZE_Z);
+  }
+
+  //number of steps in simulation
+  if (m->argc > cur_arg) {
+    finalStepCount=atoi(m->argv[cur_arg++]);
+    CkPrintf("Final Step Count:%d\n",finalStepCount);
+  }
+
+  //initializing the 3D cell array
+  cellArray = CProxy_Cell::ckNew(cellArrayDimX,cellArrayDimY,cellArrayDimZ);
+  CkPrintf("\nCells: %d X %d X %d .... created\n", cellArrayDimX, cellArrayDimY, cellArrayDimZ);
+
+  //initializing the 6D compute array
+  computeArray = CProxy_Compute::ckNew();
+  for (int x=0; x<cellArrayDimX; x++)
+    for (int y=0; y<cellArrayDimY; y++)
+      for (int z=0; z<cellArrayDimZ; z++)
+        cellArray(x, y, z).createComputes();
+
+  thisProxy.run();
+  delete m;
+}
+
+//constructor for chare object migration
+Main::Main(CkMigrateMessage* msg): CBase_Main(msg) { 
+  __sdag_init();
+}
+
+//pup routine incase the main chare moves, pack important information
+void Main::pup(PUP::er &p) {
+  CBase_Main::pup(p);
+   __sdag_pup(p);
+}
+
+#include "leanmd.def.h"
diff --git a/examples/charm++/leanmd/Main.h b/examples/charm++/leanmd/Main.h
new file mode 100644 (file)
index 0000000..23af581
--- /dev/null
@@ -0,0 +1,21 @@
+#ifndef __MAIN_H__
+#define __MAIN_H__
+
+extern /* readonly */ CProxy_Cell cellArray;
+extern /* readonly */ CProxy_Compute computeArray;
+
+extern /* readonly */ int cellArrayDimX;
+extern /* readonly */ int cellArrayDimY;
+extern /* readonly */ int cellArrayDimZ;
+
+//central controller chare
+class Main : public CBase_Main {
+  private:
+    Main_SDAG_CODE
+
+  public:
+    Main(CkArgMsg* msg);
+    Main(CkMigrateMessage* msg);
+    void pup(PUP::er &p);
+};
+#endif
diff --git a/examples/charm++/leanmd/Makefile b/examples/charm++/leanmd/Makefile
new file mode 100644 (file)
index 0000000..b73d08d
--- /dev/null
@@ -0,0 +1,30 @@
+# to be set appropiately
+#CHARMBASE      = $(HOME)/charm/net-linux-x86_64-syncft
+CHARMBASE      = $(HOME)/collegestuff/charm/net-darwin-x86_64/
+CHARMC         = $(CHARMBASE)/bin/charmc
+
+OPTS            = -O3
+
+all: leanmd
+
+leanmd: Main.o Cell.o Compute.o leanmd.decl.h
+       $(CHARMC) $(OPTS) -module DummyLB \
+       -language charm++ -o leanmd Main.o Cell.o Compute.o
+
+Main.o: Main.cc Main.h leanmd.decl.h defs.h
+       $(CHARMC) $(OPTS) -o Main.o Main.cc
+
+Cell.o: Cell.cc Cell.h leanmd.decl.h defs.h
+       $(CHARMC) $(OPTS) -o Cell.o Cell.cc
+
+leanmd.decl.h: leanmd.ci
+       $(CHARMC) leanmd.ci
+
+Compute.o: Compute.cc Compute.h leanmd.decl.h defs.h physics.h
+       $(CHARMC) $(OPTS) -o Compute.o Compute.cc
+
+test: leanmd
+       ./charmrun +p4 ./leanmd 4 4 4 101 
+
+clean:
+       rm -f *.decl.h *.def.h *.o leanmd leanmd.prj charmrun
diff --git a/examples/charm++/leanmd/README b/examples/charm++/leanmd/README
new file mode 100644 (file)
index 0000000..b29d108
--- /dev/null
@@ -0,0 +1,28 @@
+TO COMPILE-
+Do a make in this folder
+
+TO RUN-
+Either provide no arguments or provide all
+./charmrun +p<Procs> ./leanmd <dimX,dimY,dimZ,steps>
+Example - ./charmrun +p4 ./leanmd 2 2 2 1001 
+Example - ./charmrun +p4 ./leanmd 2 2 2 1001 
+
+READ ORDER
+leanmd.ci - contains the most important aspects of the program, look for run() 
+function for each chare which is called right after constructor call for each
+chare. It contains most of the parallel flow and control flow information.
+
+Cell.* - contains the implementation of important functions attached with
+Cell; focus on sendPositions() and updateProperties()
+
+Compute.* - contains the implementation for Compute; the interact() function
+is important as it does the force computation.
+
+Main.* - contains the constructor for Main chare; also the starting point 
+of the application; main chare passes control to run() of Main chare and
+constructors of Cell.
+
+physics.h - contains physics equation implementation for force computation; also
+contains the contribute for reduction of forces.
+
+defs.h - #define for constant values and some class declaration
diff --git a/examples/charm++/leanmd/defs.h b/examples/charm++/leanmd/defs.h
new file mode 100644 (file)
index 0000000..48230a4
--- /dev/null
@@ -0,0 +1,91 @@
+
+#ifndef __DEFS__
+#define __DEFS__
+
+#include "pup.h"
+
+#define HYDROGEN_MASS         (1.67 * pow( 10.0,-24))
+#define VDW_A                              (1.60694452 * pow(10.0, -134))
+#define VDW_B                              (1.031093844 * pow(10.0, -77))
+
+#define ENERGY_VAR                   (1.0 * pow(10.0,-5))
+//average of next two should be what you want as you atom density
+#define PARTICLES_PER_CELL_START  300
+#define PARTICLES_PER_CELL_END         300
+
+
+#define DEFAULT_DELTA         1        // in femtoseconds
+
+#define DEFAULT_FIRST_LDB     20
+#define DEFAULT_LDB_PERIOD    20
+#define DEFAULT_FT_PERIOD     100000
+
+#define KAWAY_X               1
+#define KAWAY_Y               1
+#define KAWAY_Z               1
+#define NBRS_X               (2*KAWAY_X+1)
+#define NBRS_Y                     (2*KAWAY_Y+1)
+#define NBRS_Z                     (2*KAWAY_Z+1)
+#define NUM_NEIGHBORS          (NBRS_X * NBRS_Y * NBRS_Z)
+
+#define CELLARRAY_DIM_X              3
+#define CELLARRAY_DIM_Y              3
+#define CELLARRAY_DIM_Z              3
+#define PTP_CUT_OFF                    12 // cut off for atom to atom interactions
+#define CELL_MARGIN                    4  // constant diff between cutoff and cell size
+#define CELL_SIZE_X                    (PTP_CUT_OFF + CELL_MARGIN)/KAWAY_X
+#define CELL_SIZE_Y                    (PTP_CUT_OFF + CELL_MARGIN)/KAWAY_Y
+#define CELL_SIZE_Z                    (PTP_CUT_OFF + CELL_MARGIN)/KAWAY_Z
+#define CELL_ORIGIN_X                0
+#define CELL_ORIGIN_Y                0
+#define CELL_ORIGIN_Z                0
+
+#define MIGRATE_STEPCOUNT          20
+#define DEFAULT_FINALSTEPCOUNT 1001
+#define MAX_VELOCITY                 30.0
+
+#define WRAP_X(a)              (((a)+cellArrayDimX)%cellArrayDimX)
+#define WRAP_Y(a)              (((a)+cellArrayDimY)%cellArrayDimY)
+#define WRAP_Z(a)              (((a)+cellArrayDimZ)%cellArrayDimZ)
+
+struct vec3 {
+  double x, y, z;
+
+  vec3(double d = 0.0) : x(d), y(d), z(d) { }
+  vec3(double x_, double y_, double z_) : x(x_), y(y_), z(z_) { }
+
+  inline vec3& operator += (const vec3 &rhs) {
+    x += rhs.x; y += rhs.y; z += rhs.z;
+    return *this;
+  }
+  inline vec3& operator -= (const vec3 &rhs) {
+    return *this += (rhs * -1.0);
+  }
+  inline vec3 operator* (const double d) const {
+    return vec3(d*x, d*y, d*z);
+  }
+  inline vec3 operator- (const vec3& rhs) const {
+    return vec3(x - rhs.x, y - rhs.y, z - rhs.z);
+  }
+};
+inline double dot(const vec3& a, const vec3& b) {
+  return a.x*b.x + a.y*b.y + a.z*b.z;
+}
+PUPbytes(vec3)
+
+//class for keeping track of the properties for a particle
+struct Particle {
+  double mass;
+  //   Position, acceleration, velocity
+  vec3 pos,acc,vel,force;
+
+  // Function for pupping properties
+  void pup(PUP::er &p) {
+    p | mass;
+    p | pos;
+    p | acc;
+    p | vel;
+       p | force;
+  }
+};
+#endif
diff --git a/examples/charm++/leanmd/leanmd.ci b/examples/charm++/leanmd/leanmd.ci
new file mode 100644 (file)
index 0000000..b41f5dc
--- /dev/null
@@ -0,0 +1,98 @@
+mainmodule leanmd {
+
+  readonly CProxy_Main mainProxy;       //central controller
+  readonly CProxy_Cell cellArray;     //array that houses atoms
+  readonly CProxy_Compute computeArray; //computational kernels
+
+  readonly int cellArrayDimX;          // X dimension of the Cell array
+  readonly int cellArrayDimY;          // Y dimension of the Cell array
+  readonly int cellArrayDimZ;          // Y dimension of the Cell array
+  readonly int finalStepCount;         // number of steps in the simulaion
+
+  //central controller chare
+  mainchare [migratable] Main {
+    entry Main(CkArgMsg* msg);
+    entry [reductiontarget] void computesCreated(); //called when computes have been created
+    entry [reductiontarget] void done(); //called when computes have been created
+
+    //this get called after the main constructor has called cell array constructor
+    entry void run() {
+      //compute array created, call for section creation
+      when computesCreated() atomic {
+        computeArray.doneInserting();
+        CkPrintf("Starting simulation .... \n\n");
+        cellArray.run();
+        computeArray.run();
+      }
+      //when all chares are done
+      when done() atomic {
+        CkPrintf("Simulation completed\n");
+        CkExit();
+      }
+    };
+  };
+
+  //message used to convey particle to computes
+  message ParticleDataMsg{
+    vec3 part[];
+  };
+
+  //chares to house atoms
+  array [3D] Cell {
+    entry Cell();  
+    entry void createComputes();    //call to insert computes that I need  
+    entry void receiveParticles(CkVec<Particle> updates);   //receive atoms that have migrated to neighboring cells to me
+    entry void receiveForces(int iter, vec3 forces[n], int n);   //receives forces from my computes on my atoms
+
+    //function to perform iterations for Cells
+    entry void run() {
+      for(stepCount = 1; stepCount <= finalStepCount; stepCount++) {
+        //send current atom positions to my computes 
+        atomic { sendPositions(); } 
+
+        //update properties of atoms using new force values 
+        // when reduceForces(vec3 forces[n], int n) atomic { updateProperties(forces, n); }
+        for(forceCount=0; forceCount < inbrs; forceCount++){
+               when receiveForces[stepCount](int iter, vec3 forces[n], int n) atomic { addForces(forces); }
+               }
+
+               atomic { updateProperties(); }
+
+        if ((stepCount %  MIGRATE_STEPCOUNT) == 0) {
+          //send atoms that have moved beyond my cell to neighbors
+          atomic { migrateParticles(); } 
+    
+          //receive particles from my neighbors
+          for(updateCount = 0; updateCount < inbrs; updateCount++) {
+            when receiveParticles(CkVec<Particle> &updates) atomic {
+              for(int i=0; i < updates.length(); i++) {
+                particles.push_back(updates[i]);    //add particles that have moved from neighboring cells to my cell
+              }
+            }
+          }
+        }
+      }
+      //everything done, reduction on kinetic energy
+      atomic { contribute(0,NULL,CkReduction::nop,CkCallback(CkReductionTarget(Main,done),mainProxy)); }
+    };
+  };
+
+  //chares that do force computations for pair of cells
+  array [6D] Compute {
+    entry Compute();
+    entry void calculateForces(ParticleDataMsg *msg);
+
+    entry void run() {
+      for(stepCount = 1; stepCount <= finalStepCount; stepCount++) {
+        //self interaction check
+        if(thisIndex.x1==thisIndex.x2 && thisIndex.y1==thisIndex.y2 && thisIndex.z1==thisIndex.z2) {
+          when calculateForces[stepCount](ParticleDataMsg *msg) atomic { selfInteract(msg); }
+        } else {
+          //receive positions from two cells - buffer the first one
+          when calculateForces[stepCount](ParticleDataMsg *msg) atomic { bufferedMsg = msg; }
+          when calculateForces[stepCount](ParticleDataMsg *msg) atomic { interact(msg); }
+        }
+      }
+    };
+  };
+};
diff --git a/examples/charm++/leanmd/physics.h b/examples/charm++/leanmd/physics.h
new file mode 100644 (file)
index 0000000..4366606
--- /dev/null
@@ -0,0 +1,116 @@
+#ifndef __PHYSICS_H__
+#define __PHYSICS_H__
+
+
+extern /* readonly */ CkGroupID mCastGrpID;
+
+extern /* readonly */ int cellArrayDimX;       // Number of Chare Rows
+extern /* readonly */ int cellArrayDimY;       // Number of Chare Columns
+extern /* readonly */ int cellArrayDimZ;
+extern /* readonly */ int finalStepCount; 
+
+#define BLOCK_SIZE     512
+
+//function to calculate forces among 2 lists of atoms
+inline void calcPairForces(ParticleDataMsg* first, ParticleDataMsg* second, int stepCount) {
+  int i, j, jpart, ptpCutOffSqd, diff;
+  int firstLen = first->lengthAll;
+  int secondLen = second->lengthAll;
+  double powTwenty, powTen, r, rsqd, f, fr;
+  vec3 separation, force;
+  double rSix, rTwelve;
+
+  vec3 *firstmsg = new vec3[firstLen];
+  vec3 *secondmsg = new vec3[secondLen];
+  //check for wrap around and adjust locations accordingly
+  if (abs(first->x - second->x) > 1){
+    diff = CELL_SIZE_X * cellArrayDimX;
+    if (second->x < first->x)
+      diff = -1 * diff; 
+    for (i = 0; i < firstLen; i++)
+      first->part[i].x += diff;
+  }
+  if (abs(first->y - second->y) > 1){
+    diff = CELL_SIZE_Y * cellArrayDimY;
+    if (second->y < first->y)
+      diff = -1 * diff; 
+    for (i = 0; i < firstLen; i++)
+      first->part[i].y += diff;
+  }
+  if (abs(first->z - second->z) > 1){
+    diff = CELL_SIZE_Z * cellArrayDimZ;
+    if (second->z < first->z)
+      diff = -1 * diff; 
+    for (i = 0; i < firstLen; i++)
+      first->part[i].z += diff;
+  } 
+  ptpCutOffSqd = PTP_CUT_OFF * PTP_CUT_OFF;
+  powTen = pow(10.0, -10);
+  powTwenty = pow(10.0, -20);
+
+  int i1, j1;
+  for(i1 = 0; i1 < firstLen; i1=i1+BLOCK_SIZE)
+    for(j1 = 0; j1 < secondLen; j1=j1+BLOCK_SIZE)
+      for(i = i1; i < i1+BLOCK_SIZE && i < firstLen; i++) {
+        for(jpart = j1; jpart < j1+BLOCK_SIZE && jpart < secondLen; jpart++) {
+          separation = first->part[i] - second->part[jpart];
+          rsqd = dot(separation, separation);
+          if (rsqd >= 0.001 && rsqd < ptpCutOffSqd) {
+            rsqd = rsqd * powTwenty;
+            r = sqrt(rsqd);
+            rSix = ((double)rsqd) * rsqd * rsqd;
+            rTwelve = rSix * rSix;
+            f = (double)(VDW_A / rTwelve - VDW_B / rSix);
+            fr = f /r;
+            force = separation * (fr * powTen);
+            firstmsg[i] += force;
+            secondmsg[jpart] -= force;
+          }
+        }
+      }
+
+  cellArray(first->x, first->y, first->z).receiveForces(stepCount,firstmsg,firstLen);
+  cellArray(second->x, second->y, second->z).receiveForces(stepCount,secondmsg,secondLen);
+
+  delete firstmsg;
+  delete secondmsg;
+  delete first;
+  delete second;
+}
+
+//function to calculate forces among atoms in a single list
+inline void calcInternalForces(ParticleDataMsg* first, int stepCount) {
+  int i, j, ptpCutOffSqd;
+  int firstLen = first->lengthAll;
+  double powTwenty, powTen, firstx, firsty, firstz, rx, ry, rz, r, rsqd, fx, fy, fz, f, fr;
+  vec3 firstpos, separation, force;
+  double rSix, rTwelve;
+  vec3 *firstmsg = new vec3[firstLen];
+
+  ptpCutOffSqd = PTP_CUT_OFF * PTP_CUT_OFF;
+  powTen = pow(10.0, -10);
+  powTwenty = pow(10.0, -20);
+  for(i = 0; i < firstLen; i++){
+    firstpos = first->part[i];
+    for(j = i+1; j < firstLen; j++) {
+      // computing base values
+      separation = firstpos - first->part[j];
+      rsqd = dot(separation, separation);
+      if(rsqd >= 0.001 && rsqd < ptpCutOffSqd){
+        rsqd = rsqd * powTwenty;
+        r = sqrt(rsqd);
+        rSix = ((double)rsqd) * rsqd * rsqd;
+        rTwelve = rSix * rSix;
+        f = (double)(VDW_A / rTwelve - VDW_B / rSix);
+        fr = f /r;
+        force = separation * (fr * powTen);
+        firstmsg[j] += force;
+        firstmsg[i] -= force;
+      }
+    }
+  }
+  cellArray(first->x, first->y, first->z).receiveForces(stepCount,firstmsg,firstLen);
+  delete firstmsg;
+  delete first;
+}
+#endif