Example charm++/matmul: simplify logic for when to pass blocks
[charm.git] / examples / charm++ / Molecular2D / Compute.C
1 /** \file Compute.h
2  *  Author: Abhinav S Bhatele
3  *  Date Created: July 1st, 2008
4  *
5  */
6
7 #include "common.h"
8 #ifdef RUN_LIVEVIZ
9   #include "liveViz.h"
10 #endif
11 #include "Patch.decl.h"
12 #include "Compute.h"
13
14 extern /* readonly */ CProxy_Main mainProxy;
15 extern /* readonly */ CProxy_Patch patchArray;
16 extern /* readonly */ CProxy_Compute computeArray;
17
18 extern /* readonly */ int numParts;
19 extern /* readonly */ int patchArrayDimX;       // Number of Chare Rows
20 extern /* readonly */ int patchArrayDimY;       // Number of Chare Columns
21 extern /* readonly */ int patchSize;
22 extern /* readonly */ double radius;
23 extern /* readonly */ int finalStepCount; 
24 extern /* readonly */ double stepTime; 
25
26 extern double A;                        // Force Calculation parameter 1
27 extern double B;                        // Force Calculation parameter 2
28
29 // Compute - Default constructor
30 Compute::Compute() {
31   cellCount = 0;
32   bufferedX = 0;
33   bufferedY = 0;
34 }
35
36 Compute::Compute(CkMigrateMessage *msg) { }
37   
38 // Function to receive vector of particles
39 void Compute::interact(CkVec<Particle> particles, int x, int y ) {
40
41   int i;
42
43   // self interaction check
44   if( thisIndex.x == thisIndex.z && thisIndex.y == thisIndex.w ) {
45     interact(particles,particles);
46     patchArray( x, y).updateForces(particles);
47   } else {
48     if(cellCount == 0) {
49       bufferedX = x;
50       bufferedY = y;
51       bufferedParticles = particles;
52       cellCount++;
53     } else if (cellCount == 1) {
54       // if both particle sets are received, compute interaction
55       cellCount = 0;
56       interact(bufferedParticles,particles);
57       patchArray(bufferedX, bufferedY).updateForces(bufferedParticles);
58       patchArray(x, y).updateForces(particles);
59     }
60   }
61 }
62
63 // Function to compute all the interactions between pairs of particles in two sets
64 void Compute::interact(CkVec<Particle> &first, CkVec<Particle> &second){
65   int i, j;
66   for(i = 0; i < first.length(); i++)
67     for(j = 0; j < second.length(); j++)
68       interact(first[i], second[j]);
69 }
70
71 // Function for computing interaction among two particles
72 // There is an extra test for interaction of identical particles, in which case there is no effect
73 void Compute::interact(Particle &first, Particle &second){
74   float rx,ry,rz,r,fx,fy,fz,f;
75
76   // computing base values
77   rx = first.x - second.x;
78   ry = first.y - second.y;
79   r = sqrt(rx*rx + ry*ry);
80
81   // We include 0.000001 to ensure that r doesn't tend to zero in the force calculation
82   // if(r < 0.000001 || r >= DEFAULT_RADIUS)
83   if(r < 0.000001 || r >= patchSize)
84     return;
85
86   f = A / pow(r,12) - B / pow(r,6);
87   fx = f * rx / r;
88   fy = f * ry / r;
89
90   // updating particle properties
91   second.fx -= fx;
92   second.fy -= fy;
93   first.fx += fx;
94   first.fy += fy;
95 }
96