-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathphysics.h
162 lines (150 loc) · 5.65 KB
/
physics.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
#ifndef __PHYSICS_H__
#define __PHYSICS_H__
#ifdef USE_SECTION_MULTICAST
#include "ckmulticast.h"
#endif
extern /* readonly */ CProxy_Main mainProxy;
extern /* readonly */ CProxy_Patch patchArray;
extern /* readonly */ CProxy_Compute computeArray;
extern /* readonly */ CkGroupID mCastGrpID;
extern /* readonly */ int patchArrayDimX; // Number of Chare Rows
extern /* readonly */ int patchArrayDimY; // Number of Chare Columns
extern /* readonly */ int patchArrayDimZ;
extern /* readonly */ int finalStepCount;
extern /* readonly */ BigReal stepTime;
#define BLOCK_SIZE 512
inline void calcPairForces(ParticleDataMsg* first, ParticleDataMsg* second, CkSectionInfo* cookie1, CkSectionInfo* cookie2) {
int i, j, jpart, ptpCutOffSqd, diff;
int firstLen = first->lengthAll;
int secondLen = second->lengthAll;
BigReal powTwenty, rx, ry, rz, r, rsqd, fx, fy, fz, f, fr, eField, constants;
double rSix, rTwelve;
ParticleForceMsg *firstmsg = new (firstLen) ParticleForceMsg;
ParticleForceMsg *secondmsg = new (secondLen) ParticleForceMsg;
firstmsg->lengthUpdates = firstLen;
secondmsg->lengthUpdates = secondLen;
//check for wrap around and adjust locations accordingly
if (abs(first->x - second->x) > 1){
diff = PATCH_SIZE_X * patchArrayDimX;
if (second->x < first->x)
diff = -1 * diff;
for (i = 0; i < firstLen; i++)
first->part[i].coord.x += diff;
}
if (abs(first->y - second->y) > 1){
diff = PATCH_SIZE_Y * patchArrayDimY;
if (second->y < first->y)
diff = -1 * diff;
for (i = 0; i < firstLen; i++)
first->part[i].coord.y += diff;
}
if (abs(first->z - second->z) > 1){
diff = PATCH_SIZE_Z * patchArrayDimZ;
if (second->z < first->z)
diff = -1 * diff;
for (i = 0; i < firstLen; i++)
first->part[i].coord.z += diff;
}
ptpCutOffSqd = PTP_CUT_OFF * PTP_CUT_OFF;
powTwenty = pow(10.0, -20);
constants = COULOMBS_CONSTANT * ELECTRON_CHARGE * ELECTRON_CHARGE;
memset(firstmsg->forces, 0, firstLen * 3*sizeof(BigReal));
memset(secondmsg->forces, 0, secondLen * 3*sizeof(BigReal));
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++) {
eField = first->part[i].charge;
for(jpart = j1; jpart < j1+BLOCK_SIZE && jpart < secondLen; jpart++) {
rx = first->part[i].coord.x - second->part[jpart].coord.x;
ry = first->part[i].coord.y - second->part[jpart].coord.y;
rz = first->part[i].coord.z - second->part[jpart].coord.z;
rsqd = rx*rx + ry*ry + rz*rz;
if (rsqd >= 0.001 && rsqd < ptpCutOffSqd){
rsqd = rsqd * powTwenty;
r = sqrt(rsqd);
rSix = ((double)rsqd) * rsqd * rsqd;
rTwelve = rSix * rSix;
f = (BigReal)(VDW_A / rTwelve - VDW_B / rSix);
f -= eField * constants * second->part[jpart].charge / rsqd;
fr = f /r;
fx = rx * fr;
fy = ry * fr;
fz = rz * fr;
secondmsg->forces[jpart].x -= fx;
secondmsg->forces[jpart].y -= fy;
secondmsg->forces[jpart].z -= fz;
firstmsg->forces[i].x += fx;
firstmsg->forces[i].y += fy;
firstmsg->forces[i].z += fz;
}
}
}
#ifdef USE_SECTION_MULTICAST
CkMulticastMgr *mCastGrp = CProxy_CkMulticastMgr(mCastGrpID).ckLocalBranch();
CkGetSectionInfo(*cookie1, first);
mCastGrp->contribute(sizeof(BigReal)*3*firstmsg->lengthUpdates, firstmsg->forces, CkReduction::sum_double, *cookie1);
CkGetSectionInfo(*cookie2, second);
mCastGrp->contribute(sizeof(BigReal)*3*secondmsg->lengthUpdates, secondmsg->forces, CkReduction::sum_double, *cookie2);
delete firstmsg;
delete secondmsg;
#else
patchArray(first->x, first->y, first->z).receiveForces(firstmsg);
patchArray(second->x, second->y, second->z).receiveForces(secondmsg);
#endif
delete first;
delete second;
}
// Local function to compute all the interactions between pairs of particles in two sets
inline void calcInternalForces(ParticleDataMsg* first, CkSectionInfo *cookie1) {
int i, j, ptpCutOffSqd;
int firstLen = first->lengthAll;
BigReal powTwenty, firstx, firsty, firstz, rx, ry, rz, r, rsqd, fx, fy, fz, f, fr, eField, constants;
double rSix, rTwelve;
ParticleForceMsg *firstmsg = new (firstLen) ParticleForceMsg;
firstmsg->lengthUpdates = firstLen;
constants = COULOMBS_CONSTANT * ELECTRON_CHARGE * ELECTRON_CHARGE;
memset(firstmsg->forces, 0, firstLen * 3*sizeof(BigReal));
ptpCutOffSqd = PTP_CUT_OFF * PTP_CUT_OFF;
powTwenty = pow(10.0, -20);
for(i = 0; i < firstLen; i++){
eField = first->part[i].charge;
firstx = first->part[i].coord.x;
firsty = first->part[i].coord.y;
firstz = first->part[i].coord.z;
for(j = i+1; j < firstLen; j++) {
// computing base values
rx = firstx - first->part[j].coord.x;
ry = firsty - first->part[j].coord.y;
rz = firstz - first->part[j].coord.z;
rsqd = rx*rx + ry*ry + rz*rz;
if(rsqd >= 0.001 && rsqd < ptpCutOffSqd){
rsqd = rsqd * powTwenty;
r = sqrt(rsqd);
rSix = ((double)rsqd) * rsqd * rsqd;
rTwelve = rSix * rSix;
f = (BigReal)(VDW_A / rTwelve - VDW_B / rSix);
f -= eField * constants * first->part[j].charge / (rsqd);
fr = f /r;
fx = rx * fr;
fy = ry * fr;
fz = rz * fr;
firstmsg->forces[j].x += fx;
firstmsg->forces[j].y += fy;
firstmsg->forces[j].z += fz;
firstmsg->forces[i].x -= fx;
firstmsg->forces[i].y -= fy;
firstmsg->forces[i].z -= fz;
}
}
}
#ifdef USE_SECTION_MULTICAST
CkMulticastMgr *mCastGrp = CProxy_CkMulticastMgr(mCastGrpID).ckLocalBranch();
mCastGrp->contribute(sizeof(BigReal)*3*firstmsg->lengthUpdates, firstmsg->forces, CkReduction::sum_double, *cookie1);
delete firstmsg;
#else
patchArray(first->x, first->y, first->z).receiveForces(firstmsg);
#endif
delete first;
}
#endif