97fa43e8a53bddc5da9c9229421281267fea31e5
[charm.git] / examples / charm++ / cell / md / patch.C
1 #include "patch.h"
2 #include "selfCompute.h"
3 #include "pairCompute.h"
4 #include "main.h"
5
6
7 void Patch::randomizeParticles() {
8
9   // Fill in a box with electrons that initially have no velocity
10   for (int i = 0; i < numParticles; i++) {
11     particleX[i] = (randf() * SIM_BOX_SIDE_LEN) - (SIM_BOX_SIDE_LEN / two);
12     particleY[i] = (randf() * SIM_BOX_SIDE_LEN) - (SIM_BOX_SIDE_LEN / two);
13     particleZ[i] = (randf() * SIM_BOX_SIDE_LEN) - (SIM_BOX_SIDE_LEN / two);
14     particleQ[i] = ELECTRON_CHARGE;
15     particleM[i] = ELECTRON_MASS;
16     velocityX[i] = zero;
17     velocityY[i] = zero;
18     velocityZ[i] = zero;
19
20     // DMK - DEBUG
21     #if DUMP_INITIAL_PARTICLE_DATA != 0
22       CkPrintf("[INFO] :: Patch[%02d,%02d,%02d]::randomizeParticles() - particle[%04d] = { "
23                "px:%+6e, py:%+6e, pz:%+6e, q:%+6e, m:%+6e, vx:%+6e, vy:%+6e, vz:%+6e }\n",
24                thisIndex.x, thisIndex.y, thisIndex.z, i,
25                particleX[i], particleY[i], particleZ[i], particleQ[i], particleM[i],
26                velocityX[i], velocityY[i], velocityZ[i]
27               );
28     #endif
29       //EJB DEBUG
30     #if SANITY_CHECK
31       CkAssert(finite(particleX[i]));
32       CkAssert(finite(particleY[i]));
33       CkAssert(finite(particleZ[i]));
34       CkAssert(finite(velocityX[i]));
35       CkAssert(finite(velocityY[i]));
36       CkAssert(finite(velocityZ[i]));
37     #endif 
38
39   }
40 }
41
42
43 MD_FLOAT Patch::randf() {
44   const int mask = 0x7FFFFFFF;
45   return (((MD_FLOAT)(rand() % mask)) / ((MD_FLOAT)(mask)));
46
47 }
48
49
50 Patch::Patch() {
51   particleX = particleY = particleZ = NULL;
52   particleQ = particleM = NULL;
53   forceSumX = forceSumY = forceSumZ = NULL;
54   velocityX = velocityY = velocityZ = NULL;
55   numParticles = -1;
56 }
57
58
59 Patch::Patch(CkMigrateMessage* msg) {
60   CkAbort("Patch::Patch(CkMigrateMessage* msg) not implemented yet");
61 }
62
63
64 Patch::~Patch() {
65   if (particleX != NULL) { CmiFreeAligned(particleX); particleX = NULL; }
66   if (particleY != NULL) { CmiFreeAligned(particleY); particleY = NULL; }
67   if (particleZ != NULL) { CmiFreeAligned(particleZ); particleZ = NULL; }
68   if (particleQ != NULL) { CmiFreeAligned(particleQ); particleQ = NULL; }
69   if (particleM != NULL) { CmiFreeAligned(particleM); particleM = NULL; }
70   if (forceSumX != NULL) { CmiFreeAligned(forceSumX); forceSumX = NULL; }
71   if (forceSumY != NULL) { CmiFreeAligned(forceSumY); forceSumY = NULL; }
72   if (forceSumZ != NULL) { CmiFreeAligned(forceSumZ); forceSumZ = NULL; }
73   if (velocityX != NULL) { CmiFreeAligned(velocityX); velocityX = NULL; }
74   if (velocityY != NULL) { CmiFreeAligned(velocityY); velocityY = NULL; }
75   if (velocityZ != NULL) { CmiFreeAligned(velocityZ); velocityZ = NULL; }
76   numParticles = 0;
77 }
78
79
80 void Patch::init(int numParticles) {
81   init_common(numParticles);
82 }
83
84 void Patch::init(int numParticles, CProxy_ProxyPatch proxyPatchProxy) {
85   #if USE_PROXY_PATCHES != 0
86     this->proxyPatchProxy = proxyPatchProxy;
87   #endif
88   init_common(numParticles);
89 }
90
91 void Patch::init_common(int numParticles) {
92
93   // Allocate memory for the particles
94   this->numParticles = numParticles;
95   particleX = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
96   particleY = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
97   particleZ = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
98   particleQ = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
99   particleM = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
100   forceSumX = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
101   forceSumY = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
102   forceSumZ = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
103   velocityX = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
104   velocityY = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
105   velocityZ = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
106
107   // Initialize the particles
108   randomizeParticles();
109
110   // Create an array section for the pair computes this patch interacts with
111   #if USE_ARRAY_SECTIONS != 0
112
113     // Enumerate the pair computes that this patch interacts with
114     const int patchIndex = PATCH_XYZ_TO_I(thisIndex.x, thisIndex.y, thisIndex.z);
115     const int numPatches = numPatchesX * numPatchesY * numPatchesZ;
116
117     // Lower section
118     CkVec<CkArrayIndex2D> pairIndexes_lower;
119     for (int i = 0; i < patchIndex; i++) {
120       pairIndexes_lower.push_back(CkArrayIndex2D(i, patchIndex));
121     }
122     pairComputeArraySection_lower = CProxySection_PairCompute::ckNew(pairComputeArrayProxy, pairIndexes_lower.getVec(), pairIndexes_lower.size());
123
124     // Upper section
125     CkVec<CkArrayIndex2D> pairIndexes_upper;
126     for (int i = patchIndex + 1; i < numPatches; i++) {
127       pairIndexes_upper.push_back(CkArrayIndex2D(patchIndex, i));
128     }
129     pairComputeArraySection_upper = CProxySection_PairCompute::ckNew(pairComputeArrayProxy, pairIndexes_upper.getVec(), pairIndexes_upper.size());
130
131   #endif
132
133   // Check in with the main proxy
134   mainProxy.initCheckIn();
135 }
136
137
138 void Patch::startIteration() { startIteration_common(1); }
139 void Patch::startIterations(int numIters) { startIteration_common(numIters); }
140 void Patch::startIteration_common(int numIters) {
141
142   // Set the number of remaining time stops
143   CkAssert(numIters > 0);
144   remainingIterations = numIters;
145
146   // Reset the number of expected computes that will check in
147   //   NOTE: pair compute for every other patch, self compute for this patch
148   remainingForceCheckIns = numPatchesX * numPatchesY * numPatchesZ;
149
150   // Clear the force sum arrays
151   #if 1
152     memset(forceSumX, 0, sizeof(MD_FLOAT) * numParticles);
153     memset(forceSumY, 0, sizeof(MD_FLOAT) * numParticles);
154     memset(forceSumZ, 0, sizeof(MD_FLOAT) * numParticles);
155   #else
156     register MD_VEC* fsx = (MD_VEC*)forceSumX;
157     register MD_VEC* fsy = (MD_VEC*)forceSumY;
158     register MD_VEC* fsz = (MD_VEC*)forceSumZ;
159     const MD_VEC zero_vec = vspread_MDF(zero);
160     register const int numParticles_vec = numParticles / myvec_numElems;
161     for (int i = 0; i < numParticles_vec; i++) {
162       fsx[i] = zero_vec;
163       fsy[i] = zero_vec;
164       fsz[i] = zero_vec;
165     }
166   #endif
167
168   // DMK - DEBUG
169   NetworkProgress
170
171   // Send particle data to pair computes
172
173   #if USE_PROXY_PATCHES != 0
174
175     proxyPatchProxy.patchData(numParticles, particleX, particleY, particleZ, particleQ);
176
177   #elif USE_ARRAY_SECTIONS != 0
178
179     // Send particle data to self computes
180     selfComputeArrayProxy(thisIndex.x, thisIndex.y, thisIndex.z).patchData(numParticles, particleX, particleY, particleZ, particleQ);
181
182     pairComputeArraySection_lower.patchData(numParticles, particleX, particleY, particleZ, particleQ, 1);
183     pairComputeArraySection_upper.patchData(numParticles, particleX, particleY, particleZ, particleQ, 0);
184
185   #else
186
187     // Send particle data to self computes
188     selfComputeArrayProxy(thisIndex.x, thisIndex.y, thisIndex.z).patchData(numParticles, particleX, particleY, particleZ, particleQ);
189
190     const int index = PATCH_XYZ_TO_I(thisIndex.x, thisIndex.y, thisIndex.z);
191     for (int i = 0; i < index; i++) {
192
193       // DMK - DEBUG
194       NetworkProgress
195
196       pairComputeArrayProxy(i, index).patchData(numParticles, particleX, particleY, particleZ, particleQ, 1);
197     }
198     const int numPatches = numPatchesX * numPatchesY * numPatchesZ;
199     for (int i = index + 1; i < numPatches; i++) {
200
201       // DMK - DEBUG
202       NetworkProgress
203
204       pairComputeArrayProxy(index, i).patchData(numParticles, particleX, particleY, particleZ, particleQ, 0);
205     }
206
207     // DMK - DEBUG
208     NetworkProgress
209
210   #endif
211 }
212
213
214 void Patch::forceCheckIn(int numParticles, MD_FLOAT* forceX, MD_FLOAT* forceY, MD_FLOAT* forceZ) {
215   forceCheckIn(numParticles, forceX, forceY, forceZ, 1);
216 }
217 void Patch::forceCheckIn(int numParticles, MD_FLOAT* forceX, MD_FLOAT* forceY, MD_FLOAT* forceZ, int numForceCheckIns) {
218
219   // Accumulate the force data
220   #if 0
221     register MD_VEC* fsx = (MD_VEC*)forceSumX;
222     register MD_VEC* fsy = (MD_VEC*)forceSumY;
223     register MD_VEC* fsz = (MD_VEC*)forceSumZ;
224     register MD_VEC* fx = (MD_VEC*)forceX;
225     register MD_VEC* fy = (MD_VEC*)forceY;
226     register MD_VEC* fz = (MD_VEC*)forceZ;
227     register const int numParticles_vec = numParticles / myvec_numElems;
228     register int i;
229     for (i = 0; i < numParticles_vec; i++) {
230       fsx[i] = vadd_MDF(fsx[i], fx[i]);
231       fsy[i] = vadd_MDF(fsy[i], fy[i]);
232       fsz[i] = vadd_MDF(fsz[i], fz[i]);
233     }
234   #else
235     for (int i = 0; i < numParticles; i++) {
236 #if SANITY_CHECK
237       CkAssert(finite(forceX[i]));
238       CkAssert(finite(forceY[i]));
239       CkAssert(finite(forceZ[i]));
240 #endif 
241       forceSumX[i] += forceX[i];
242       forceSumY[i] += forceY[i];
243       forceSumZ[i] += forceZ[i];
244     }
245   #endif
246
247   // Count the incoming forced data and integrate if all force data has arrived
248   remainingForceCheckIns -= numForceCheckIns;
249   if (remainingForceCheckIns <= 0) {
250     thisProxy(thisIndex.x, thisIndex.y, thisIndex.z).integrate();
251   }
252 }
253
254
255 void Patch::integrate_callback() {
256
257   // DMK - DEBUG
258   #if ENABLE_USER_EVENTS != 0
259     double __start_time__ = CmiWallTimer();
260   #endif
261
262   // DMK - DEBUG
263   #if COUNT_FLOPS != 0
264     globalFlopCount += localFlopCount;
265   #endif
266
267   // DMK - DEBUG
268   NetworkProgress
269
270   // Decrement the counter containing the number of remaining iterations.  If there are
271   //   more iterations, do another one, otherwise, check in with main
272   remainingIterations--;
273   if (remainingIterations > 0) {
274     startIteration_common(remainingIterations);
275   } else {
276     mainProxy.patchCheckIn();
277   }
278
279   // DMK - DEBUG
280   NetworkProgress
281
282   // DMK - DEBUG
283   #if ENABLE_USER_EVENTS != 0
284     double __end_time__ = CmiWallTimer();
285     traceUserBracketEvent(PROJ_USER_EVENT_PATCH_INTEGRATE_CALLBACK, __start_time__, __end_time__);
286   #endif
287 }
288
289
290 ProxyPatch::ProxyPatch(int patchIndex) {
291   this->patchIndex = patchIndex;
292   numParticles = -1;
293   particleX = particleY = particleZ = particleQ = NULL;
294   forceSumX = forceSumY = forceSumZ = NULL;
295 }
296
297
298 ProxyPatch::ProxyPatch(CkMigrateMessage *msg) {
299   CkAbort("ProxyPatch::ProxyPatch(CkMigrateMessage *msg) not implemented yet");
300 }
301
302
303 ProxyPatch::~ProxyPatch() {
304   if (particleX != NULL) { CmiFreeAligned(particleX); particleX = NULL; }
305   if (particleY != NULL) { CmiFreeAligned(particleY); particleY = NULL; }
306   if (particleZ != NULL) { CmiFreeAligned(particleZ); particleZ = NULL; }
307   if (particleQ != NULL) { CmiFreeAligned(particleQ); particleQ = NULL; }
308   if (forceSumX != NULL) { CmiFreeAligned(forceSumX); forceSumX = NULL; }
309   if (forceSumY != NULL) { CmiFreeAligned(forceSumY); forceSumY = NULL; }
310   if (forceSumZ != NULL) { CmiFreeAligned(forceSumZ); forceSumZ = NULL; }
311   numParticles = -1;
312 }
313
314
315 void ProxyPatch::init(int numParticles) {
316
317   // Allocate memory for the particles
318   this->numParticles = numParticles;
319   particleX = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
320   particleY = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
321   particleZ = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
322   particleQ = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
323   forceSumX = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
324   forceSumY = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
325   forceSumZ = (MD_FLOAT*)(CmiMallocAligned(numParticles * sizeof(MD_FLOAT), 128));
326
327   // Check in with the main proxy
328   mainProxy.initCheckIn();
329 }
330
331
332 void ProxyPatch::patchData(int numParticles, MD_FLOAT* particleX, MD_FLOAT* particleY, MD_FLOAT* particleZ, MD_FLOAT* particleQ) {
333
334   // Copy in the updated particle data
335   memcpy(this->particleX, particleX, numParticles * sizeof(MD_FLOAT));
336   memcpy(this->particleY, particleY, numParticles * sizeof(MD_FLOAT));
337   memcpy(this->particleZ, particleZ, numParticles * sizeof(MD_FLOAT));
338   memcpy(this->particleQ, particleQ, numParticles * sizeof(MD_FLOAT));
339
340   // Clear out the force arrays
341   memset(this->forceSumX, 0, numParticles * sizeof(MD_FLOAT));
342   memset(this->forceSumY, 0, numParticles * sizeof(MD_FLOAT));
343   memset(this->forceSumZ, 0, numParticles * sizeof(MD_FLOAT));
344
345   // Reset the patch checkin counters
346   checkInCount = 0;
347
348   // Call patchData on the local self compute
349   const int patchX = PATCH_I_TO_X(patchIndex);
350   const int patchY = PATCH_I_TO_Y(patchIndex);
351   const int patchZ = PATCH_I_TO_Z(patchIndex);
352   SelfCompute* localSelfCompute = selfComputeArrayProxy(patchX, patchY, patchZ).ckLocal();
353   if (localSelfCompute != NULL) {
354     localSelfCompute->patchData(numParticles, this->particleX, this->particleY, this->particleZ, this->particleQ, thisProxy);
355     checkInCount++;
356   }
357
358   // Call patchData on all local pair computes
359   const int myPe = CkMyPe();
360   for (int i = 0; i < patchIndex; i++) {
361     PairCompute* localPairCompute = pairComputeArrayProxy(i, patchIndex).ckLocal();
362     if (localPairCompute != NULL) {
363       localPairCompute->patchData(numParticles, this->particleX, this->particleY, this->particleZ, this->particleQ, 1, thisProxy);
364       checkInCount++;
365     }
366   }
367   const int numPatches = numPatchesX * numPatchesY * numPatchesZ;
368   for (int i = patchIndex + 1; i < numPatches; i++) {
369     PairCompute* localPairCompute = pairComputeArrayProxy(patchIndex, i).ckLocal();
370     if (localPairCompute != NULL) {
371       localPairCompute->patchData(numParticles, this->particleX, this->particleY, this->particleZ, this->particleQ, 0, thisProxy);
372       checkInCount++;
373     }
374   }
375
376   numToCheckIn = checkInCount;
377 }
378
379
380 void ProxyPatch::forceCheckIn(int numParticles, MD_FLOAT* forceX, MD_FLOAT* forceY, MD_FLOAT* forceZ) {
381
382   // Accumulate the force data
383   #if USE_PROXY_PATCHES != 0  // Calls will be local and pointers will be aligned, so take advantage and vectorize the code
384     register MD_VEC* forceX_vec = (MD_VEC*)forceX;
385     register MD_VEC* forceY_vec = (MD_VEC*)forceY;
386     register MD_VEC* forceZ_vec = (MD_VEC*)forceZ;
387     register MD_VEC* forceSumX_vec = (MD_VEC*)forceSumX;
388     register MD_VEC* forceSumY_vec = (MD_VEC*)forceSumY;
389     register MD_VEC* forceSumZ_vec = (MD_VEC*)forceSumZ;
390     const int numParticles_vec = numParticles / myvec_numElems;
391     for (int i = 0; i < numParticles_vec; i++) {
392       forceSumX_vec[i] += forceX_vec[i];
393       forceSumY_vec[i] += forceY_vec[i];
394       forceSumZ_vec[i] += forceZ_vec[i];
395     }
396   #else
397     for (int i = 0; i < numParticles; i++) {
398       forceSumX[i] += forceX[i];
399       forceSumY[i] += forceY[i];
400       forceSumZ[i] += forceZ[i];
401     }
402   #endif
403
404   // Once all computes this proxy called have contributed forces, send the data back to the patch itself
405   checkInCount--;
406   if (checkInCount <= 0) {
407     const int patchX = PATCH_I_TO_X(patchIndex);
408     const int patchY = PATCH_I_TO_Y(patchIndex);
409     const int patchZ = PATCH_I_TO_Z(patchIndex);
410     patchArrayProxy(patchX, patchY, patchZ).forceCheckIn(numParticles, forceSumX, forceSumY, forceSumZ, numToCheckIn);
411   }
412 }
413
414
415 #include "patch.def.h"