Merge branch 'charm' into amicableErrorMessages
[charm.git] / examples / multiphaseSharedArrays / moldyn / moldyn.C
1 // -*- mode: c++; tab-width: 4 -*-
2
3 // When running 1D, make NEPP = COL1
4 // When running 2D, same
5 // When running 3D, make NEPP = subset of COL1
6
7 #include "nepp.h"
8 #include "msa/msa.h"
9
10 class XYZ { // coords, forces
11 public:
12     double x, y, z;
13     XYZ() { x = y = z = 0.0; }
14     XYZ(const int rhs) { x = y = z = (double)rhs; } // identity value
15     XYZ& operator+= (const XYZ& rhs) { x += rhs.x; y += rhs.y; z += rhs.z; return *this;}
16     XYZ& negate() { x = -x; y=-y; z=-z; return *this;}
17 };
18 PUPbytes(XYZ);
19
20 class AtomInfo {
21 public:
22     double mass, charge;
23   AtomInfo() : mass(0.0), charge(0.0) { }
24     AtomInfo(const int rhs) { mass = charge = (double)rhs; } // identity value
25     AtomInfo& operator+= (const AtomInfo& rhs) { } // we're not calling accumulate on this
26 };
27 PUPbytes(AtomInfo);
28
29 typedef MSA::MSA1D<XYZ, DefaultEntry<XYZ,false>, NEPP> XyzMSA;
30 typedef MSA::MSA1D<AtomInfo, DefaultEntry<AtomInfo,false>, NEPP> AtomInfoMSA;
31 typedef MSA::MSA2D<bool, DefaultEntry<bool,false>, NEPP, MSA_ROW_MAJOR> NeighborMSA;
32
33 #include "moldyn.decl.h"
34
35 #include <assert.h>
36 #include <math.h>
37 #include "params.h"
38
39 const double epsilon = 0.00000001;
40 inline int notequal(double v1, double v2)
41 {
42     return (fabs(v1 - v2) > epsilon);
43 }
44
45 class moldyn : public CBase_moldyn
46 {
47 protected:
48     double start_time;
49     CProxy_WorkerArray workers;
50     int reallyDone;
51
52 public:
53     moldyn(CkArgMsg* m)
54     {
55         // Usage: a.out [number_of_worker_threads [max_bytes]]
56         if(m->argc >1 ) NUM_WORKERS=atoi(m->argv[1]);
57         if(m->argc >2 ) NUM_ATOMS=atoi(m->argv[2]);
58         if(m->argc >3 ) CACHE_SIZE_BYTES=atoi(m->argv[3]);
59         if(m->argc >4 ) detailedTimings= ((atoi(m->argv[4])!=0)?true:false) ; // 1D, 2D, 3D
60         delete m;
61         reallyDone = 0;
62
63         XyzMSA coords(NUM_ATOMS, NUM_WORKERS, CACHE_SIZE_BYTES);
64         XyzMSA forces(NUM_ATOMS, NUM_WORKERS, CACHE_SIZE_BYTES);
65         AtomInfoMSA atominfo(NUM_ATOMS, NUM_WORKERS, CACHE_SIZE_BYTES);
66         NeighborMSA nbrList(NUM_ATOMS, NUM_ATOMS, NUM_WORKERS, CACHE_SIZE_BYTES);
67
68         workers = CProxy_WorkerArray::ckNew(coords, forces, atominfo, nbrList, NUM_WORKERS, NUM_WORKERS);
69         workers.ckSetReductionClient(new CkCallback(CkIndex_moldyn::done(NULL), thisProxy));
70
71         start_time = CkWallTimer();
72         workers.Start();
73     }
74
75     // This method gets called twice, and should only terminate the
76     // second time.
77     void done(CkReductionMsg* m)
78     {
79         int *ip = (int*)m->getData();
80         bool prefetchWorked = (*ip==0);
81         delete m;
82
83         if (reallyDone == 0) {
84             workers.Kontinue();
85             reallyDone++;
86
87             double end_time = CkWallTimer();
88
89             const char TAB = '\t';
90
91             char hostname[100];
92             gethostname(hostname, 100);
93
94             ckout << CkNumPes() << TAB
95                   << NUM_WORKERS << TAB
96                   << "nepp " << NEPP << TAB
97                                   << "atom " << NUM_ATOMS << TAB
98                   << end_time - start_time << TAB
99                   << CACHE_SIZE_BYTES << TAB
100                   << (runPrefetchVersion? (prefetchWorked?"Y":"N"): "U") << " "
101                   << hostname
102                   << endl;
103
104         } else {
105             CkExit();
106         }
107     }
108 };
109
110 // Returns start and end
111 void GetMyIndices(unsigned int maxIndex, unsigned int myNum, unsigned int numWorkers,
112                   unsigned int& start, unsigned int& end)
113 {
114     int rangeSize = maxIndex / numWorkers;
115     if(myNum < maxIndex % numWorkers)
116     {
117         start = myNum * (rangeSize + 1);
118         end = start + rangeSize;
119     }
120     else
121     {
122         start = myNum * rangeSize + maxIndex % numWorkers;
123         end = start + rangeSize - 1;
124     }
125 }
126
127 class WorkerArray : public CBase_WorkerArray
128 {
129 private:
130     // prefetchWorked keeps track of whether the prefetches succeeded or not.
131     bool prefetchWorked;
132     CkVec<double> times;
133     CkVec<const char*> description;
134
135     // ================================================================
136     // 2D calculations
137
138     inline int numWorkers2D() {
139         static int n = 0;
140
141         if (n==0) {
142             n = (int)(sqrt(numWorkers));
143             CkAssert(n*n == numWorkers);
144         }
145
146         return n;
147     }
148
149     // Convert a 1D ChareArray index into a 2D x dimension index
150     inline unsigned int toX() {
151         return thisIndex/numWorkers2D();
152     }
153     // Convert a 1D ChareArray index into a 2D y dimension index
154     inline unsigned int toY() {
155         return thisIndex%numWorkers2D();
156     }
157
158     // ================================================================
159
160 protected:
161     XyzMSA coords;
162     XyzMSA forces;
163     AtomInfoMSA atominfo;
164         AtomInfoMSA::Read rAtominfo;
165     NeighborMSA nbrList;
166
167     unsigned int numAtoms, numWorkers;
168
169     void EnrollArrays()
170     {
171         coords.enroll(numWorkers);
172         forces.enroll(numWorkers);
173         atominfo.enroll(numWorkers);
174         nbrList.enroll(numWorkers);
175     }
176
177     void FillArrays()
178     {
179         /*
180         // fill in our portion of the array
181         unsigned int rowStart, rowEnd, colStart, colEnd;
182         GetMyIndices(rows1, thisIndex, numWorkers, rowStart, rowEnd);
183         GetMyIndices(cols2, thisIndex, numWorkers, colStart, colEnd);
184
185         // fill them in with 1
186         for(unsigned int r = rowStart; r <= rowEnd; r++)
187             for(unsigned int c = 0; c < cols1; c++)
188                 arr1.set(r, c) = 1.0;
189
190         for(unsigned int c = colStart; c <= colEnd; c++)
191             for(unsigned int r = 0; r < rows2; r++)
192                 arr2.set(r, c) = 1.0;
193         */
194     }
195
196     XYZ calculateForce(const XYZ &coordsi, const AtomInfo &atominfoi, const XYZ &coordsj, const AtomInfo &atominfoj)
197     {
198         XYZ result;
199         return result;
200     }
201
202     XYZ integrate(const AtomInfo &atominfok, const XYZ &forcesk)
203     {
204         XYZ result;
205         return result;
206     }
207
208     double distance(unsigned int i, unsigned int j)
209     {
210         return 0;
211     }
212
213   double distance(const XYZ &a, const XYZ &b)
214   {
215         double dx = a.x - b.x,
216           dy = a.y - b.y,
217           dz = a.z - b.z;
218         return sqrt(dx*dx + dy*dy + dz*dz);
219   }
220
221   void PlimptonMD(XyzMSA::Handle hCoords, XyzMSA::Handle hForces, NeighborMSA::Handle hNbr)
222   {
223         unsigned int i_start, i_end, j_start, j_end;
224         GetMyIndices(NUM_ATOMS-1, toX(), numWorkers2D(), i_start, i_end);
225         GetMyIndices(NUM_ATOMS-1, toY(), numWorkers2D(), j_start, j_end);
226
227         XyzMSA::Read rCoords = hCoords.syncToRead();
228         NeighborMSA::Read rNbr = hNbr.syncToRead();
229
230         for (unsigned int timestep = 0; timestep < NUM_TIMESTEPS; timestep++) {
231           // Force calculation for a section of the interaction matrix
232           XyzMSA::Accum aForces = hForces.syncToAccum();
233           for (unsigned int i = i_start; i< i_end; i++)
234                 for (unsigned int j = j_start; j< j_end; j++)
235                   if (rNbr(i,j)) {
236                         XYZ force = calculateForce(rCoords(i),
237                                                                            rAtominfo(i),
238                                                                            rCoords(j),
239                                                                            rAtominfo(j));
240                         aForces(i) += force;
241                         aForces(j) += force.negate();
242                   }
243
244           // Movement Integration for our subset of atoms
245           unsigned int myAtomsBegin, myAtomsEnd;
246           XyzMSA::Read rForces = aForces.syncToRead();
247           XyzMSA::Write wCoords = rCoords.syncToWrite();
248           for (unsigned int k = myAtomsBegin; k<myAtomsEnd; k++)
249                 wCoords(k) = integrate(rAtominfo(k), rForces(k));
250
251           // Neighbor list recalculation for our section of the interaction matrix
252           rCoords = wCoords.syncToRead();
253           if  (timestep % 8 == 0) { // update neighbor list every 8 steps
254                 NeighborMSA::Write wNbr = rNbr.syncToWrite();
255                 for (unsigned int i = i_start; i< i_end; i++)
256                   for (unsigned int j = j_start; j< j_end; j++)
257                         if (distance(rCoords(i), rCoords(j)) < CUTOFF_DISTANCE) {
258                           wNbr.set(i,j) = true;
259                           wNbr.set(j,i) = true;
260                         } else {
261                           wNbr.set(i,j) = false;
262                           wNbr.set(j,i) = false;
263                         }
264                 rNbr = wNbr.syncToRead();
265           }
266
267           hForces = rForces;
268         }
269   }
270
271     void TestResults(bool prod_test=true)
272     {
273         /*
274         int errors = 0;
275         bool ok=true;
276
277         // verify the results, print out first error only
278         ok=true;
279         for(unsigned int r = 0; ok && r < rows1; r++) {
280             for(unsigned int c = 0; ok && c < cols1; c++) {
281                 if(notequal(arr1.get(r, c), 1.0)) {
282                     ckout << "[" << CkMyPe() << "," << thisIndex << "] arr1 -- Illegal element at (" << r << "," << c << ") " << arr1.get(r,c) << endl;
283                     ok=false;
284                     errors++;
285                 }
286             }
287         }
288
289         ok=true;
290         for(unsigned int c = 0; ok && c < cols2; c++) {
291             for(unsigned int r = 0; ok && r < rows2; r++) {
292                 if(notequal(arr2.get(r, c), 1.0)) {
293                     ckout << "[" << CkMyPe() << "," << thisIndex << "] arr2 -- Illegal element at (" << r << "," << c << ") " << arr2.get(r,c) << endl;
294                     ok=false;
295                     errors++;
296                 }
297             }
298         }
299
300         //arr1.FreeMem();
301         //arr2.FreeMem();
302
303         if(prod_test)
304         {
305             ok = true;
306             for(unsigned int c = 0; ok && c < cols2; c++) {
307                 for(unsigned int r = 0; ok && r < rows1; r++) {
308                     if(notequal(prod.get(r,c), 1.0 * cols1)) {
309                         ckout << "[" << CkMyPe() << "] result  -- Illegal element at (" << r << "," << c << ") " << prod.get(r,c) << endl;
310                         ok=false;
311                         errors++;
312                     }
313                 }
314             }
315         }
316
317         if (errors!=0) CkAbort("Incorrect array elements detected!");
318         */
319     }
320
321     void Contribute()
322     {
323         int dummy = prefetchWorked?0:1;
324         contribute(sizeof(int), &dummy, CkReduction::sum_int);
325     }
326
327 public:
328     WorkerArray(const XyzMSA &coords_, const XyzMSA &forces_, AtomInfoMSA &atominfo_,
329                 NeighborMSA &nbrList_, unsigned int numWorkers_)
330         : coords(coords_), forces(forces_), atominfo(atominfo_), nbrList(nbrList_),
331           numWorkers(numWorkers_), prefetchWorked(false), numAtoms(coords.length())
332     {
333         // ckout << "w" << thisIndex << ":" << rows1 << " " << cols1 << " " << cols2 << endl;
334         times.push_back(CkWallTimer());
335         description.push_back("constr");
336     }
337
338     WorkerArray(CkMigrateMessage* m) {}
339
340     ~WorkerArray()
341     {
342     }
343
344     void Start()
345     {
346         times.push_back(CkWallTimer()); // 1
347         description.push_back("   start");
348
349         EnrollArrays();
350         times.push_back(CkWallTimer()); // 2
351         description.push_back("   enroll");
352
353         if(verbose) ckout << thisIndex << ": filling" << endl;
354         FillArrays();
355         times.push_back(CkWallTimer()); // 3
356         description.push_back("  fill");
357
358         if(verbose) ckout << thisIndex << ": syncing" << endl;
359         //SyncArrays();
360         times.push_back(CkWallTimer()); // 4
361         description.push_back("    sync");
362
363         if (do_test) TestResults(0);
364
365         if(verbose) ckout << thisIndex << ": product" << endl;
366
367         PlimptonMD(coords.getInitialWrite(), forces.getInitialWrite(), nbrList.getInitialWrite());
368         times.push_back(CkWallTimer()); // 5
369         description.push_back("    work");
370
371         Contribute();
372     }
373
374     void Kontinue()
375     {
376         times.push_back(CkWallTimer()); // 6
377         description.push_back("    redn");
378
379         if(verbose) ckout << thisIndex << ": testing" << endl;
380         if (do_test) TestResults();
381         times.push_back(CkWallTimer()); // 7
382         description.push_back("    test");
383         Contribute();
384
385         if (detailedTimings) {
386             if (thisIndex == 0) {
387                 for(int i=1; i<description.length(); i++)
388                     ckout << description[i] << " ";
389                 ckout << endl;
390             }
391             ckout << "w" << thisIndex << ":";
392             for(int i=1; i<times.length(); i++)
393                 ckout << times[i]-times[i-1] << " ";
394             ckout << endl;
395         }
396     }
397 };
398
399 #include "moldyn.def.h"