1/*
2 * vector_dist_performance_common.hpp
3 *
4 * Created on: Dec 25, 2016
5 * Author: i-bird
6 */
7
8#ifndef SRC_VECTOR_PERFORMANCE_VECTOR_DIST_PERFORMANCE_COMMON_HPP_
9#define SRC_VECTOR_PERFORMANCE_VECTOR_DIST_PERFORMANCE_COMMON_HPP_
10
11#include "Vector/vector_dist.hpp"
12
13/*! \brief Calculate and put particles' forces
14 *
15 * \param NN Cell list
16 * \param vd Distributed vector
17 * \param r_cut Cut-off radius
18 */
19template<unsigned int dim, size_t prp = 0, typename T, typename V> void calc_forces(T & NN, V & vd, float r_cut)
20{
21 auto it_v = vd.getDomainIterator();
22
23 float sum[dim];
24
25 for (size_t i = 0; i < dim; i++)
26 sum[i] = 0;
27
28 while (it_v.isNext())
29 {
30 //key
31 vect_dist_key_dx key = it_v.get();
32
33 // Get the position of the particles
34 Point<dim,float> p = vd.getPos(key);
35
36 for (size_t i = 0; i < dim; i++)
37 sum[i] = 0;
38
39 // Get the neighborhood of the particle
40 auto cell_it = NN.template getNNIterator<NO_CHECK>(NN.getCell(p));
41
42 while(cell_it.isNext())
43 {
44 auto nnp = cell_it.get();
45
46 // p != q
47 if (nnp == key.getKey())
48 {
49 ++cell_it;
50 continue;
51 }
52
53 Point<dim,float> q = vd.getPos(nnp);
54
55 if (p.distance2(q) < r_cut*r_cut)
56 {
57 //Calculate the forces
58 float num[dim];
59 for (size_t i = 0; i < dim; i++)
60 num[i] = vd.getPos(key)[i] - vd.getPos(nnp)[i];
61
62 float denom = 0;
63 for (size_t i = 0; i < dim; i++)
64 denom += num[i] * num[i];
65
66 float res[dim];
67 for (size_t i = 0; i < dim; i++)
68 res[i] = num[i] / denom;
69
70 for (size_t i = 0; i < dim; i++)
71 sum[i] += res[i];
72 }
73 //Next particle in a cell
74 ++cell_it;
75 }
76
77 //Put the forces
78 for (size_t i = 0; i < dim; i++)
79 vd.template getProp<prp>(key)[i] += sum[i];
80
81 //Next particle in cell list
82 ++it_v;
83 }
84}
85
86/*! \brief For each particle of vd calculate the accumulation of the distances of the neighborhood
87 * particles inside vd2
88 *
89 *
90 * \param NN Cell list vd
91 * \param NN2 Cell list vd2
92 * \param vd Distributed vector
93 * \param vd2 Distributed vector 2
94 * \param r_cut Cut-off radius
95 *
96 */
97template<unsigned int dim, unsigned int prp, typename T, typename V> void cross_calc(T & NN, T & NN2, V & vd, V & vd2)
98{
99 auto it_v = vd.getDomainIterator();
100
101 while (it_v.isNext())
102 {
103 //key
104 vect_dist_key_dx key = it_v.get();
105
106 // Get the position of the particles
107 Point<dim,float> p = vd.getPos(key);
108
109 // Get the neighborhood of the particle
110 auto cell_it = NN2.template getNNIterator<NO_CHECK>(NN2.getCell(p));
111
112 double sum = 0.0;
113
114 while(cell_it.isNext())
115 {
116 auto nnp = cell_it.get();
117
118 Point<dim,float> q = vd2.getPos(nnp);
119
120 sum += norm(p - q);
121
122 //Next particle in a cell
123 ++cell_it;
124 }
125
126 vd.template getProp<prp>(key) = sum;
127
128 //Next particle in cell list
129 ++it_v;
130 }
131}
132
133/*! \brief Initialize a distributed vector
134 *
135 * \param vd Distributed vector
136 */
137template<unsigned int dim, typename v_dist>
138void vd_initialize(v_dist & vd, Vcluster<HeapMemory> & v_cl)
139{
140 // The random generator engine
141 std::default_random_engine eg(v_cl.getProcessUnitID()*4313);
142 std::uniform_real_distribution<float> ud(0.0f, 1.0f);
143
144 //! [Create a vector of random elements on each processor 2D]
145
146 auto it = vd.getIterator();
147
148 while (it.isNext())
149 {
150 auto key = it.get();
151
152 for (size_t i = 0; i < dim; i++)
153 {vd.getPos(key)[i] = ud(eg);}
154
155 ++it;
156 }
157
158 vd.map();
159}
160
161/*! \brief Initialize a distributed vector
162 *
163 * \param vd Distributed vector
164 * \param box where to initialize
165 * \param npart number of particles to add
166 *
167 */
168template<unsigned int dim, typename v_dist>
169void vd_initialize_box_nomap(v_dist & vd, const Box<dim,typename v_dist::stype> & box, Vcluster<HeapMemory> & v_cl,int start, int stop)
170{
171 // The random generator engine
172 std::default_random_engine eg(v_cl.getProcessUnitID()*4313);
173 std::uniform_real_distribution<float> ud(0.0f, 1.0f);
174
175 auto it = vd.getIterator(start,stop);
176
177 while (it.isNext())
178 {
179 auto key = it.get();
180
181 for (size_t i = 0; i < dim; i++)
182 {vd.getPos(key)[i] = (box.getHigh(i) - box.getLow(i))*ud(eg) + box.getLow(i);}
183
184 ++it;
185 }
186}
187
188
189/*! \brief Initialize 2 distributed vectors with equally positioned particles
190 *
191 * \param vd, vd2 Distributed vectors
192 * \param v_cl Global vcluster
193 * \param k_int Number of particles
194 */
195template<unsigned int dim, typename v_dist> void vd_initialize_double(v_dist & vd,v_dist & vd2, Vcluster<> & v_cl, size_t k_int)
196{
197 // The random generator engine
198 std::default_random_engine eg(v_cl.getProcessUnitID()*4313);
199 std::uniform_real_distribution<float> ud(0.0f, 1.0f);
200
201 //! [Create a vector of random elements on each processor 2D]
202
203 auto it = vd.getIterator();
204
205 while (it.isNext())
206 {
207 auto key = it.get();
208
209 for (size_t i = 0; i < dim; i++)
210 {
211 vd.getPos(key)[i] = ud(eg);
212 vd2.getPos(key)[i] = vd.getPos(key)[i];
213 }
214
215 ++it;
216 }
217
218 vd.map();
219 vd2.map();
220}
221
222
223
224/*! \brief Calculate and put particles' forces
225 *
226 * \param NN Cell list hilbert
227 * \param vd Distributed vector
228 * \param r_cut Cut-off radius
229 */
230template<unsigned int dim, size_t prp = 0, typename T, typename V> void calc_forces_hilb(T & NN, V & vd, float r_cut)
231{
232 auto it_cl = NN.getIterator();
233
234 float sum[dim];
235
236 for (size_t i = 0; i < dim; i++)
237 sum[i] = 0;
238
239 while (it_cl.isNext())
240 {
241 //key
242 auto key = it_cl.get();
243
244 // Get the position of the particles
245 Point<dim,float> p = vd.getPos(key);
246
247 for (size_t i = 0; i < dim; i++)
248 sum[i] = 0;
249
250 // Get the neighborhood of the particle
251 auto cell_it = NN.template getNNIterator<NO_CHECK>(NN.getCell(p));
252
253 while(cell_it.isNext())
254 {
255 auto nnp = cell_it.get();
256
257 // p != q
258 if (nnp == key)
259 {
260 ++cell_it;
261 continue;
262 }
263
264 Point<dim,float> q = vd.getPos(nnp);
265
266 if (p.distance2(q) < r_cut*r_cut)
267 {
268 //Calculate the forces
269 float num[dim];
270 for (size_t i = 0; i < dim; i++)
271 num[i] = vd.getPos(key)[i] - vd.getPos(nnp)[i];
272
273 float denom = 0;
274 for (size_t i = 0; i < dim; i++)
275 denom += num[i] * num[i];
276
277 float res[dim];
278 for (size_t i = 0; i < dim; i++)
279 res[i] = num[i] / denom;
280
281 for (size_t i = 0; i < dim; i++)
282 sum[i] += res[i];
283 }
284 //Next particle in a cell
285 ++cell_it;
286 }
287
288 //Put the forces
289 for (size_t i = 0; i < dim; i++)
290 vd.template getProp<prp>(key)[i] += sum[i];
291
292 //Next particle in cell list
293 ++it_cl;
294 }
295}
296
297#endif /* SRC_VECTOR_PERFORMANCE_VECTOR_DIST_PERFORMANCE_COMMON_HPP_ */
298