Trace for micropic✔
  1. Preprocessing loop contracts✔
  2. make local copy of field✔
    1. Matrix.local_name_tile ~var:"fieldAtCorners" ~elem_ty:vect ~uninit_post:true ~mark_load:"loadField" ~local_var:"lFieldAtCorners" [ctx; cFor "idStep"];✔
  3. inline helper functions and reveal record fields✔
    1. Function.inline_multi [ctx; cFuns ["cornerInterpolationCoeff"; "matrix_vect_mul"; "vect_add"; "vect_mul"]];✔
    2. Record.split_fields ~typs:[particle; vect] [tSpanSeq [ctx]];✔
    3. Record.to_variables [ctx; cVarDefs ["fieldAtPos"; "pos2"; "speed2"; "accel"]];✔
  4. scale field and particles✔
    1. Variable.insert ~name:"fieldFactor" ~value:(trm_mul (trm_mul deltaT deltaT) (trm_div (find_var "pCharge") (find_var "pMass"))) [ctx; tBefore; cVarDef "lFieldAtCorners"]; (* tFirst *)✔
    2. List.iter scaleFieldAtPos dims;✔
    3. List.iter scaleSpeed2 dims;✔
    4. List.iter scaleFieldAtCorners dims;✔
    5. List.iter scaleParticles dims;✔
    6. List.iter Loop.fusion_targets [[cMark "partsPrep"]; [cMark "partsPostp"]];✔
  5. inline variables and simplify arithmetic✔
    1. Variable.unfold ~at:[cFor "idStep"] [cVarDef "fieldFactor"];✔
    2. Variable.inline [ctx; cVarDefs (Tools.cart_prod (^) ["accel_"; "pos2_"] dims)];✔
    3. Arith.(simpls_rec [expand; gather_rec]) [ctx];✔
  6. final polish✔
    1. Loop.hoist_alloc ~indep:["idStep"; "idPart"] ~dest:[tBefore; cFor "idStep"] [cVarDef "coeffs"];✔
    2. Cleanup.std (); )✔
tmp/{before8d8f79.cpp → afterc54b72.cpp} RENAMED
@@ -1,116 +1,100 @@
1
  #include <optitrust.h>
2
 
3
  typedef struct {
4
  double x;
5
  double y;
6
  double z;
7
  } vect;
8
 
9
  template <typename T, typename U>
10
  U* __struct_access_x(T* v) {
11
- __pure();
12
  __admitted();
13
  }
14
 
15
  template <typename T, typename U>
16
  U __struct_get_x(T v) {
17
- __pure();
18
  __admitted();
19
  }
20
 
21
  template <typename T, typename U>
22
  U* __struct_access_y(T* v) {
23
- __pure();
24
  __admitted();
25
  }
26
 
27
  template <typename T, typename U>
28
  U __struct_get_y(T v) {
29
- __pure();
30
  __admitted();
31
  }
32
 
33
  template <typename T, typename U>
34
  U* __struct_access_z(T* v) {
35
- __pure();
36
  __admitted();
37
  }
38
 
39
  template <typename T, typename U>
40
  U __struct_get_z(T v) {
41
- __pure();
42
  __admitted();
43
  }
44
 
45
  vect vect_add(vect v1, vect v2) {
46
- __pure();
47
  __admitted();
48
  return (vect){v1.x + v2.x, v1.y + v2.y, v1.z + v2.z};
49
  }
50
 
51
  vect vect_mul(double d, vect v) {
52
- __pure();
53
  __admitted();
54
  return (vect){d * v.x, d * v.y, d * v.z};
55
  }
56
 
57
  typedef struct {
58
  vect pos;
59
  vect speed;
60
  } particle;
61
 
62
  template <typename T, typename U>
63
  U* __struct_access_pos(T* v) {
64
- __pure();
65
  __admitted();
66
  }
67
 
68
  template <typename T, typename U>
69
  U __struct_get_pos(T v) {
70
- __pure();
71
  __admitted();
72
  }
73
 
74
  template <typename T, typename U>
75
  U* __struct_access_speed(T* v) {
76
- __pure();
77
  __admitted();
78
  }
79
 
80
  template <typename T, typename U>
81
  U __struct_get_speed(T v) {
82
- __pure();
83
  __admitted();
84
  }
85
 
86
  template <typename T, typename U>
87
  U* __struct_access_charge(T* v) {
88
- __pure();
89
  __admitted();
90
  }
91
 
92
  template <typename T, typename U>
93
  U __struct_get_charge(T v) {
94
- __pure();
95
  __admitted();
96
  }
97
 
98
  template <typename T, typename U>
99
  U* __struct_access_mass(T* v) {
100
- __pure();
101
  __admitted();
102
  }
103
 
104
  template <typename T, typename U>
105
  U __struct_get_mass(T v) {
106
- __pure();
107
  __admitted();
108
  }
109
 
110
  const double areaX = 10.;
111
 
112
  const double areaY = 10.;
113
 
114
  const double areaZ = 10.;
115
 
116
  const int gridX = 64;
@@ -121,140 +105,119 @@ const int gridZ = 64;
121
 
122
  const int nbCells = gridX * gridY * gridZ;
123
 
124
  const double cellX = areaX / gridX;
125
 
126
  const double cellY = areaY / gridY;
127
 
128
  const double cellZ = areaZ / gridZ;
129
 
130
  int int_of_double(double a) {
131
- __pure();
132
  __admitted();
133
  return (int)a - (a < 0.);
134
  }
135
 
136
  double relativePosX(double x) {
137
- __pure();
138
  __admitted();
139
  int iX = int_of_double(x / cellX);
140
  return (x - iX * cellX) / cellX;
141
  }
142
 
143
  double relativePosY(double y) {
144
- __pure();
145
  __admitted();
146
  int iY = int_of_double(y / cellY);
147
  return (y - iY * cellY) / cellY;
148
  }
149
 
150
  double relativePosZ(double z) {
151
- __pure();
152
  __admitted();
153
  int iZ = int_of_double(z / cellZ);
154
  return (z - iZ * cellZ) / cellZ;
155
  }
156
 
157
  const int nbCorners = 8;
158
 
159
  void cornerInterpolationCoeff(vect pos, double* r) {
160
- __writes("r ~> Matrix1(nbCorners)");
161
  const double rX = relativePosX(pos.x);
162
  const double rY = relativePosY(pos.y);
163
  const double rZ = relativePosZ(pos.z);
164
  const double cX = 1. + -1. * rX;
165
  const double cY = 1. + -1. * rY;
166
  const double cZ = 1. + -1. * rZ;
167
- __ghost(
168
- [&]() {
169
- __consumes("_Uninit(r ~> Matrix1(nbCorners))");
170
- __produces("_Uninit(&r[MINDEX1(8, 0)] ~> Cell)");
171
- __produces("_Uninit(&r[MINDEX1(8, 1)] ~> Cell)");
172
- __produces("_Uninit(&r[MINDEX1(8, 2)] ~> Cell)");
173
- __produces("_Uninit(&r[MINDEX1(8, 3)] ~> Cell)");
174
- __produces("_Uninit(&r[MINDEX1(8, 4)] ~> Cell)");
175
- __produces("_Uninit(&r[MINDEX1(8, 5)] ~> Cell)");
176
- __produces("_Uninit(&r[MINDEX1(8, 6)] ~> Cell)");
177
- __produces("_Uninit(&r[MINDEX1(8, 7)] ~> Cell)");
178
- __admitted();
179
- },
180
- "");
181
- r[MINDEX1(8, 0)] = cX * cY * cZ;
182
- r[MINDEX1(8, 1)] = cX * cY * rZ;
183
- r[MINDEX1(8, 2)] = cX * rY * cZ;
184
- r[MINDEX1(8, 3)] = cX * rY * rZ;
185
- r[MINDEX1(8, 4)] = rX * cY * cZ;
186
- r[MINDEX1(8, 5)] = rX * cY * rZ;
187
- r[MINDEX1(8, 6)] = rX * rY * cZ;
188
- r[MINDEX1(8, 7)] = rX * rY * rZ;
189
- __ghost(
190
- [&]() {
191
- __consumes("&r[MINDEX1(8, 0)] ~> Cell");
192
- __consumes("&r[MINDEX1(8, 1)] ~> Cell");
193
- __consumes("&r[MINDEX1(8, 2)] ~> Cell");
194
- __consumes("&r[MINDEX1(8, 3)] ~> Cell");
195
- __consumes("&r[MINDEX1(8, 4)] ~> Cell");
196
- __consumes("&r[MINDEX1(8, 5)] ~> Cell");
197
- __consumes("&r[MINDEX1(8, 6)] ~> Cell");
198
- __consumes("&r[MINDEX1(8, 7)] ~> Cell");
199
- __produces("r ~> Matrix1(nbCorners)");
200
- __admitted();
201
- },
202
- "");
203
  }
204
 
205
  vect matrix_vect_mul(double* coeffs, vect* matrix) {
206
- __reads("coeffs ~> Matrix1(nbCorners)");
207
- __reads("matrix ~> Matrix1(nbCorners)");
208
  vect res = {0., 0., 0.};
209
  for (int k = 0; k < nbCorners; k++) {
210
- __xreads("&coeffs[MINDEX1(nbCorners, k)] ~> Cell");
211
- __xreads("&matrix[MINDEX1(nbCorners, k)] ~> Cell");
212
- res = vect_add(res, vect_mul(coeffs[MINDEX1(nbCorners, k)],
213
- matrix[MINDEX1(nbCorners, k)]));
214
  }
215
  __admitted();
216
  return res;
217
  }
218
 
219
  void simulate_single_cell(double deltaT, particle* particles, int nbParticles,
220
  vect* fieldAtCorners, int nbSteps, double pCharge,
221
  double pMass) {
222
- __modifies("particles ~> Matrix1(nbParticles)");
223
- __reads("fieldAtCorners ~> Matrix1(nbCorners)");
 
 
 
 
 
 
 
 
 
 
 
224
  for (int idStep = 0; idStep < nbSteps; idStep++) {
225
  for (int idPart = 0; idPart < nbParticles; idPart++) {
226
- __xmodifies("&particles[MINDEX1(nbParticles, idPart)] ~> Cell");
227
- __ghost(
228
- [&]() {
229
- __consumes("&particles[MINDEX1(nbParticles, idPart)] ~> Cell");
230
- __produces("&particles[MINDEX1(nbParticles, idPart)].pos ~> Cell");
231
- __produces(
232
- "&particles[MINDEX1(nbParticles, idPart)].speed ~> Cell");
233
- __admitted();
234
- },
235
- "");
236
- double* const coeffs = (double*)MALLOC1(nbCorners, sizeof(double));
237
- cornerInterpolationCoeff(particles[MINDEX1(nbParticles, idPart)].pos,
238
- coeffs);
239
- const vect fieldAtPos = matrix_vect_mul(coeffs, fieldAtCorners);
240
- MFREE1(nbCorners, coeffs);
241
- const vect accel = vect_mul(pCharge / pMass, fieldAtPos);
242
- const vect speed2 =
243
- vect_add(particles[MINDEX1(nbParticles, idPart)].speed,
244
- vect_mul(deltaT, accel));
245
- const vect pos2 = vect_add(particles[MINDEX1(nbParticles, idPart)].pos,
246
- vect_mul(deltaT, speed2));
247
- particles[MINDEX1(nbParticles, idPart)].pos = pos2;
248
- particles[MINDEX1(nbParticles, idPart)].speed = speed2;
249
- __ghost(
250
- [&]() {
251
- __consumes("&particles[MINDEX1(nbParticles, idPart)].pos ~> Cell");
 
 
 
 
 
 
 
 
 
 
 
252
- __consumes(
 
 
 
253
- "&particles[MINDEX1(nbParticles, idPart)].speed ~> Cell");
254
- __produces("&particles[MINDEX1(nbParticles, idPart)] ~> Cell");
255
- __admitted();
256
- },
257
- "");
258
  }
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
259
  }
 
260
  }
1
  #include <optitrust.h>
2
 
3
  typedef struct {
4
  double x;
5
  double y;
6
  double z;
7
  } vect;
8
 
9
  template <typename T, typename U>
10
  U* __struct_access_x(T* v) {
 
11
  __admitted();
12
  }
13
 
14
  template <typename T, typename U>
15
  U __struct_get_x(T v) {
 
16
  __admitted();
17
  }
18
 
19
  template <typename T, typename U>
20
  U* __struct_access_y(T* v) {
 
21
  __admitted();
22
  }
23
 
24
  template <typename T, typename U>
25
  U __struct_get_y(T v) {
 
26
  __admitted();
27
  }
28
 
29
  template <typename T, typename U>
30
  U* __struct_access_z(T* v) {
 
31
  __admitted();
32
  }
33
 
34
  template <typename T, typename U>
35
  U __struct_get_z(T v) {
 
36
  __admitted();
37
  }
38
 
39
  vect vect_add(vect v1, vect v2) {
 
40
  __admitted();
41
  return (vect){v1.x + v2.x, v1.y + v2.y, v1.z + v2.z};
42
  }
43
 
44
  vect vect_mul(double d, vect v) {
 
45
  __admitted();
46
  return (vect){d * v.x, d * v.y, d * v.z};
47
  }
48
 
49
  typedef struct {
50
  vect pos;
51
  vect speed;
52
  } particle;
53
 
54
  template <typename T, typename U>
55
  U* __struct_access_pos(T* v) {
 
56
  __admitted();
57
  }
58
 
59
  template <typename T, typename U>
60
  U __struct_get_pos(T v) {
 
61
  __admitted();
62
  }
63
 
64
  template <typename T, typename U>
65
  U* __struct_access_speed(T* v) {
 
66
  __admitted();
67
  }
68
 
69
  template <typename T, typename U>
70
  U __struct_get_speed(T v) {
 
71
  __admitted();
72
  }
73
 
74
  template <typename T, typename U>
75
  U* __struct_access_charge(T* v) {
 
76
  __admitted();
77
  }
78
 
79
  template <typename T, typename U>
80
  U __struct_get_charge(T v) {
 
81
  __admitted();
82
  }
83
 
84
  template <typename T, typename U>
85
  U* __struct_access_mass(T* v) {
 
86
  __admitted();
87
  }
88
 
89
  template <typename T, typename U>
90
  U __struct_get_mass(T v) {
 
91
  __admitted();
92
  }
93
 
94
  const double areaX = 10.;
95
 
96
  const double areaY = 10.;
97
 
98
  const double areaZ = 10.;
99
 
100
  const int gridX = 64;
105
 
106
  const int nbCells = gridX * gridY * gridZ;
107
 
108
  const double cellX = areaX / gridX;
109
 
110
  const double cellY = areaY / gridY;
111
 
112
  const double cellZ = areaZ / gridZ;
113
 
114
  int int_of_double(double a) {
 
115
  __admitted();
116
  return (int)a - (a < 0.);
117
  }
118
 
119
  double relativePosX(double x) {
 
120
  __admitted();
121
  int iX = int_of_double(x / cellX);
122
  return (x - iX * cellX) / cellX;
123
  }
124
 
125
  double relativePosY(double y) {
 
126
  __admitted();
127
  int iY = int_of_double(y / cellY);
128
  return (y - iY * cellY) / cellY;
129
  }
130
 
131
  double relativePosZ(double z) {
 
132
  __admitted();
133
  int iZ = int_of_double(z / cellZ);
134
  return (z - iZ * cellZ) / cellZ;
135
  }
136
 
137
  const int nbCorners = 8;
138
 
139
  void cornerInterpolationCoeff(vect pos, double* r) {
 
140
  const double rX = relativePosX(pos.x);
141
  const double rY = relativePosY(pos.y);
142
  const double rZ = relativePosZ(pos.z);
143
  const double cX = 1. + -1. * rX;
144
  const double cY = 1. + -1. * rY;
145
  const double cZ = 1. + -1. * rZ;
 
 
 
 
 
 
 
 
 
 
 
 
 
 
146
+ r[0] = cX * cY * cZ;
147
+ r[1] = cX * cY * rZ;
148
+ r[2] = cX * rY * cZ;
149
+ r[3] = cX * rY * rZ;
150
+ r[4] = rX * cY * cZ;
151
+ r[5] = rX * cY * rZ;
152
+ r[6] = rX * rY * cZ;
153
+ r[7] = rX * rY * rZ;
 
 
 
 
 
 
 
 
 
 
 
 
 
 
154
  }
155
 
156
  vect matrix_vect_mul(double* coeffs, vect* matrix) {
 
 
157
  vect res = {0., 0., 0.};
158
  for (int k = 0; k < nbCorners; k++) {
 
 
159
+ res = vect_add(res, vect_mul(coeffs[k], matrix[k]));
 
160
  }
161
  __admitted();
162
  return res;
163
  }
164
 
165
  void simulate_single_cell(double deltaT, particle* particles, int nbParticles,
166
  vect* fieldAtCorners, int nbSteps, double pCharge,
167
  double pMass) {
168
+ const int fieldFactor = deltaT * deltaT * pCharge / pMass;
169
+ vect* const lFieldAtCorners = (vect*)MALLOC1(nbCorners, sizeof(vect));
170
+ for (int i1 = 0; i1 < nbCorners; i1++) {
171
+ lFieldAtCorners[i1].x = fieldAtCorners[i1].x * fieldFactor;
172
+ lFieldAtCorners[i1].y = fieldAtCorners[i1].y * fieldFactor;
173
+ lFieldAtCorners[i1].z = fieldAtCorners[i1].z * fieldFactor;
174
+ }
175
+ for (int i1 = 0; i1 < nbParticles; i1 += 1) {
176
+ particles[i1].speed.x *= deltaT;
177
+ particles[i1].speed.y *= deltaT;
178
+ particles[i1].speed.z *= deltaT;
179
+ }
180
+ double* const coeffs = (double*)MALLOC1(nbCorners, sizeof(double));
181
  for (int idStep = 0; idStep < nbSteps; idStep++) {
182
  for (int idPart = 0; idPart < nbParticles; idPart++) {
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
183
+ const double rX = relativePosX(particles[idPart].pos.x);
 
184
+ const double rY = relativePosY(particles[idPart].pos.y);
 
 
 
185
+ const double rZ = relativePosZ(particles[idPart].pos.z);
186
+ const double cX = 1. + -1. * rX;
187
+ const double cY = 1. + -1. * rY;
188
+ const double cZ = 1. + -1. * rZ;
189
+ coeffs[0] = cX * cY * cZ;
190
+ coeffs[1] = cX * cY * rZ;
191
+ coeffs[2] = cX * rY * cZ;
192
+ coeffs[3] = cX * rY * rZ;
193
+ coeffs[4] = rX * cY * cZ;
194
+ coeffs[5] = rX * cY * rZ;
195
+ coeffs[6] = rX * rY * cZ;
196
+ coeffs[7] = rX * rY * rZ;
197
+ double fieldAtPos_x = 0.;
198
+ double fieldAtPos_y = 0.;
199
+ double fieldAtPos_z = 0.;
200
+ for (int k = 0; k < nbCorners; k++) {
201
+ fieldAtPos_x += coeffs[k] * lFieldAtCorners[k].x;
202
+ fieldAtPos_y += coeffs[k] * lFieldAtCorners[k].y;
203
+ fieldAtPos_z += coeffs[k] * lFieldAtCorners[k].z;
 
 
204
  }
205
+ const double speed2_x = particles[idPart].speed.x + fieldAtPos_x;
206
+ const double speed2_y = particles[idPart].speed.y + fieldAtPos_y;
207
+ const double speed2_z = particles[idPart].speed.z + fieldAtPos_z;
208
+ particles[idPart].pos.x += speed2_x;
209
+ particles[idPart].pos.y += speed2_y;
210
+ particles[idPart].pos.z += speed2_z;
211
+ particles[idPart].speed.x = speed2_x;
212
+ particles[idPart].speed.y = speed2_y;
213
+ particles[idPart].speed.z = speed2_z;
214
+ }
215
+ }
216
+ MFREE1(nbCorners, coeffs);
217
+ for (int i1 = 0; i1 < nbParticles; i1 += 1) {
218
+ particles[i1].speed.z /= deltaT;
219
+ particles[i1].speed.y /= deltaT;
220
+ particles[i1].speed.x /= deltaT;
221
  }
222
+ MFREE1(nbCorners, lFieldAtCorners);
223
  }
 Show: