OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
main.cpp
1 #include "Grid/grid_dist_id.hpp"
2 #include "data_type/aggregate.hpp"
3 #include "timer.hpp"
4 
9 constexpr int U = 0;
10 constexpr int V = 1;
11 constexpr int phi = 2;
12 constexpr int normal = 3;
13 constexpr int tgrad_u = 4;
14 constexpr int tgrad_v = 5;
15 constexpr int U_next = 6;
16 constexpr int V_next = 7;
17 
18 constexpr int x = 0;
19 constexpr int y = 1;
20 constexpr int z = 2;
21 
22 typedef sgrid_dist_soa<3,double,aggregate<double,double,double,double[3],double[3],double[3],double,double> > sgrid_type;
23 
24 void init(sgrid_type & grid, Box<3,double> & domain)
25 {
27 
28  auto it = grid.getGridIterator();
29  Point<3,double> p1({0.5,0.5,0.5});
30 
31  double sx = grid.spacing(0);
32 
33  Sphere<3,double> sph1(p1,0.3);
34  Sphere<3,double> sph2(p1,0.3 - sx*10);
35  Sphere<3,double> sph_zero(p1,0.3 - sx*5);
36 
37  while (it.isNext())
38  {
39  // Get the local grid key
40  auto key = it.get_dist();
41  auto keyg = it.get();
42 
43  Point<3,double> pc;
44  Point<3,double> vp;
45 
46  for (int i = 0 ; i < 3 ; i++)
47  {pc.get(i) = keyg.get(i) * it.getSpacing(i);}
48 
49  // Check if the point is in the first sphere
50  if (sph1.isInside(pc) == true && sph2.isInside(pc) == false)
51  {
52  Point<3,double> pn = pc - p1;
53  pn /= pn.norm();
54  double theta = acos(pn * Point<3,double>({0.0,0.0,1.0}));
55  Point<3,double> pn_ = pn;
56  pn_[2] = 0.0;
57  pn_ /= pn_.norm();
58  double aphi = acos(pn_ * Point<3,double>({1.0,0.0,0.0}));
59 
60  // Create a perturbation in the solid angle
61  if (theta > 0.6 && theta < 0.8 && aphi > 0.0 && aphi < 0.2)
62  {
63  grid.template insert<U>(key) = 0.5;
64  grid.template insert<V>(key) = 0.25;
65  }
66  else
67  {
68  grid.template insert<U>(key) = 1.0;
69  grid.template insert<V>(key) = 0.0;
70  }
71  grid.template insert<phi>(key) = sph_zero.distance(pc);
72  grid.template insert<normal>(key)[0] = pn[0];
73  grid.template insert<normal>(key)[1] = pn[1];
74  grid.template insert<normal>(key)[2] = pn[2];
75 
76  // Old values U and V
77  grid.template insert<U_next>(key) = 0.0;
78  grid.template insert<V_next>(key) = 0.0;
79  }
80 
81  ++it;
82  }
83 
85 }
86 
87 template<unsigned int U_src,unsigned int V_src,unsigned int U_dst, unsigned int V_dst>
88 void extend(sgrid_type & grid, size_t (& sz)[3],double (& spacing)[3])
89 {
90  double delta = 1e-10;
91  double max = 0.0;
92 
93  auto func_extend = [delta,&spacing](auto & grid, auto & ids,
94  unsigned char * mask_sum)
95  {
96  Vc::double_v phi_c;
97  Vc::double_v s;
98 
99  Vc::double_v Uext = 0.0;
100  Vc::double_v Vext = 0.0;
101 
102  Vc::double_v n[3];
103  Vc::double_v dir;
104 
105  Vc::double_v Uc;
106  Vc::double_v Vc;
107  Vc::double_v Uc_xm;
108  Vc::double_v Vc_xm;
109  Vc::double_v Uc_ym;
110  Vc::double_v Vc_ym;
111  Vc::double_v Uc_zm;
112  Vc::double_v Vc_zm;
113 
114  Vc::double_v Uc_xp;
115  Vc::double_v Vc_xp;
116  Vc::double_v Uc_yp;
117  Vc::double_v Vc_yp;
118  Vc::double_v Uc_zp;
119  Vc::double_v Vc_zp;
120 
121  load_crs<x,0,phi>(phi_c,grid,ids);
122  load_crs_v<x,0,x,normal>(n[x],grid,ids);
123  load_crs_v<x,0,y,normal>(n[y],grid,ids);
124  load_crs_v<x,0,z,normal>(n[z],grid,ids);
125 
126  load_crs<x,0,U_src>(Uc,grid,ids);
127  load_crs<x,0,V_src>(Vc,grid,ids);
128  load_crs<x,-1,U_src>(Uc_xm,grid,ids);
129  load_crs<x,-1,V_src>(Vc_xm,grid,ids);
130  load_crs<y,-1,U_src>(Uc_ym,grid,ids);
131  load_crs<y,-1,V_src>(Vc_ym,grid,ids);
132  load_crs<z,-1,U_src>(Uc_zm,grid,ids);
133  load_crs<z,-1,V_src>(Vc_zm,grid,ids);
134  load_crs<x,1,U_src>(Uc_xp,grid,ids);
135  load_crs<x,1,V_src>(Vc_xp,grid,ids);
136  load_crs<y,1,U_src>(Uc_yp,grid,ids);
137  load_crs<y,1,V_src>(Vc_yp,grid,ids);
138  load_crs<z,1,U_src>(Uc_zp,grid,ids);
139  load_crs<z,1,V_src>(Vc_zp,grid,ids);
140 
141  s = phi_c / sqrt(phi_c*phi_c + delta*delta);
142 
143  dir = s*n[0];
144  auto dir_pos = dir > 0;
145  auto dir_neg = dir < 0;
146 
147  Uext += Vc::iif(dir_pos,dir * (Uc - Uc_xm)/spacing[0],Vc::double_v(0.0));
148  Vext += Vc::iif(dir_pos,dir * (Vc - Vc_xm)/spacing[0],Vc::double_v(0.0));
149  Uext += Vc::iif(dir_neg,dir * (Uc_xp - Uc)/spacing[0],Vc::double_v(0.0));
150  Vext += Vc::iif(dir_neg,dir * (Vc_xp - Vc)/spacing[0],Vc::double_v(0.0));
151 
152  dir = s*n[1];
153  dir_pos = dir > 0;
154  dir_neg = dir < 0;
155 
156  Uext += Vc::iif(dir_pos,dir * (Uc - Uc_ym)/spacing[1],Vc::double_v(0.0));
157  Vext += Vc::iif(dir_pos,dir * (Vc - Vc_ym)/spacing[1],Vc::double_v(0.0));
158  Uext += Vc::iif(dir_neg,dir * (Uc_yp - Uc)/spacing[1],Vc::double_v(0.0));
159  Vext += Vc::iif(dir_neg,dir * (Vc_yp - Vc)/spacing[1],Vc::double_v(0.0));
160 
161  dir = s*n[2];
162  dir_pos = dir > 0;
163  dir_neg = dir < 0;
164 
165  Uext += Vc::iif(dir_pos,dir * (Uc - Uc_zm)/spacing[2],Vc::double_v(0.0));
166  Vext += Vc::iif(dir_pos,dir * (Vc - Vc_zm)/spacing[2],Vc::double_v(0.0));
167  Uext += Vc::iif(dir_neg,dir * (Uc_zp - Uc)/spacing[2],Vc::double_v(0.0));
168  Vext += Vc::iif(dir_neg,dir * (Vc_zp - Vc)/spacing[2],Vc::double_v(0.0));
169 
170  Uext = Uc - 0.0003*Uext;
171  Vext = Vc - 0.0003*Vext;
172 
173  store_crs<U_dst>(grid,Uext,ids);
174  store_crs<V_dst>(grid,Vext,ids);
175  };
176 
177  grid.template conv_cross_ids<1,double>({0,0,0},{sz[0] - 1, sz[1] - 1, sz[2] - 1},func_extend);
178 }
179 
180 int main(int argc, char* argv[])
181 {
182  openfpm_init(&argc,&argv);
183 
184  // domain
185  Box<3,double> domain({0.0,0.0,0.0},{2.5,2.5,2.5});
186 
187  // grid size
188  size_t sz[3] = {512,512,512};
189 
190  // Define periodicity of the grid
191  periodicity<3> bc = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
192 
193  // Ghost in grid unit
194  Ghost<3,long int> g(1);
195 
196  // deltaT
197  double deltaT = 0.3;
198 
199  // Diffusion constant for specie U
200  double du = 1*1e-5;
201 
202  // Diffusion constant for specie V
203  double dv = 0.5*1e-5;
204 
205 #ifdef TEST_RUN
206  // Number of timesteps
207  size_t timeSteps = 200;
208 #else
209  // Number of timesteps
210  size_t timeSteps = 100000;
211 #endif
212 
213  // K and F (Physical constant in the equation)
214  double K = 0.053;
215  double F = 0.014;
216 
217  sgrid_type grid(sz,domain,g,bc);
218 
219 
220  // spacing of the grid on x and y
221  double spacing[3] = {grid.spacing(0),grid.spacing(1),grid.spacing(2)};
222 
223  init(grid,domain);
224 
225  // sync the ghost
226  size_t count = 0;
227  grid.template ghost_get<U,V>();
228 
229  // because we assume that spacing[x] == spacing[y] we use formula 2
230  // and we calculate the prefactor of Eq 2
231  double uFactor = deltaT * du;
232  double vFactor = deltaT * dv;
233 
234  auto & v_cl = create_vcluster();
235 
236  timer tot_sim;
237  tot_sim.start();
238 
239  for (size_t i = 0; i < timeSteps ; ++i)
240  {
241  auto func_grad = [&spacing](auto & grid, auto & ids,
242  unsigned char * mask_sum){
243 
244  Vc::double_v n[3];
245 
246  Vc::double_v Uc;
247  Vc::double_v xmU;
248  Vc::double_v ymU;
249  Vc::double_v zmU;
250 
251  Vc::double_v Vc;
252  Vc::double_v xmV;
253  Vc::double_v ymV;
254  Vc::double_v zmV;
255 
256  Vc::double_v u_out[3];
257  Vc::double_v v_out[3];
258 
259  load_crs<x,-1,U>(xmU,grid,ids);
260  load_crs<y,-1,U>(ymU,grid,ids);
261  load_crs<z,-1,U>(zmU,grid,ids);
262  load_crs<x,0,U>(Uc,grid,ids);
263 
264  load_crs<x,-1,V>(xmV,grid,ids);
265  load_crs<y,-1,V>(ymV,grid,ids);
266  load_crs<z,-1,V>(zmV,grid,ids);
267  load_crs<x,0,V>(Vc,grid,ids);
268 
269  load_crs_v<x,0,x,normal>(n[x],grid,ids);
270  load_crs_v<x,0,y,normal>(n[y],grid,ids);
271  load_crs_v<x,0,z,normal>(n[z],grid,ids);
272 
273  u_out[0] = (1.0-n[0]*n[0])*(Uc - xmU)/spacing[0] + (-n[1]*n[1])*(Uc - ymU)/spacing[1] + (-n[2]*n[2])*(Uc - zmU)/spacing[2];
274  u_out[1] = (-n[0]*n[0])*(Uc - xmU)/spacing[0] + (1.0-n[1]*n[1])*(Uc - ymU)/spacing[1] + (-n[2]*n[2])*(Uc - zmU)/spacing[2];
275  u_out[2] = (-n[0]*n[0])*(Uc - xmU)/spacing[0] + (-n[1]*n[1])*(Uc - ymU)/spacing[1] + (1.0-n[2]*n[2])*(Uc - zmU)/spacing[2];
276 
277  v_out[0] = (1.0-n[0]*n[0])*(Vc - xmV)/spacing[0] + (-n[1]*n[1])*(Vc - ymV)/spacing[1] + (-n[2]*n[2])*(Vc - zmV)/spacing[2];
278  v_out[1] = (-n[0]*n[0])*(Vc - xmV)/spacing[0] + (1.0-n[1]*n[1])*(Vc - ymV)/spacing[1] + (-n[2]*n[2])*(Vc - zmV)/spacing[2];
279  v_out[2] = (-n[0]*n[0])*(Vc - xmV)/spacing[0] + (-n[1]*n[1])*(Vc - ymV)/spacing[1] + (1.0-n[2]*n[2])*(Vc - zmV)/spacing[2];
280 
281  Vc::Mask<double> surround;
282 
283  for (int i = 0 ; i < Vc::double_v::Size ; i++)
284  {surround[i] = (mask_sum[i] == 6);}
285 
286  u_out[0] = Vc::iif(surround,u_out[0],Vc::double_v(0.0));
287  u_out[1] = Vc::iif(surround,u_out[1],Vc::double_v(0.0));
288  u_out[2] = Vc::iif(surround,u_out[2],Vc::double_v(0.0));
289 
290  v_out[0] = Vc::iif(surround,v_out[0],Vc::double_v(0.0));
291  v_out[1] = Vc::iif(surround,v_out[1],Vc::double_v(0.0));
292  v_out[2] = Vc::iif(surround,v_out[2],Vc::double_v(0.0));
293 
294  store_crs_v<tgrad_u,x>(grid,u_out[0],ids);
295  store_crs_v<tgrad_u,y>(grid,u_out[1],ids);
296  store_crs_v<tgrad_u,z>(grid,u_out[2],ids);
297 
298  store_crs_v<tgrad_v,x>(grid,v_out[0],ids);
299  store_crs_v<tgrad_v,y>(grid,v_out[1],ids);
300  store_crs_v<tgrad_v,z>(grid,v_out[2],ids);
301  };
302 
303  grid.template conv_cross_ids<1,double>({0,0,0},{sz[0]-1,sz[1] - 1,sz[2] - 1},func_grad);
304 
305  auto func_lap = [&spacing,uFactor,vFactor,deltaT,K,F](auto & grid, auto & ids,
306  unsigned char * mask_sum){
307 
308  Vc::double_v gradU_px;
309  Vc::double_v gradU_py;
310  Vc::double_v gradU_pz;
311 
312  Vc::double_v gradU_x;
313  Vc::double_v gradU_y;
314  Vc::double_v gradU_z;
315 
316  Vc::double_v gradV_px;
317  Vc::double_v gradV_py;
318  Vc::double_v gradV_pz;
319 
320  Vc::double_v gradV_x;
321  Vc::double_v gradV_y;
322  Vc::double_v gradV_z;
323 
324  Vc::double_v lapU;
325  Vc::double_v lapV;
326 
327  Vc::double_v Uc;
328  Vc::double_v Vc;
329 
330  Vc::double_v outU;
331  Vc::double_v outV;
332 
333  load_crs_v<x,1,x,tgrad_u>(gradU_px,grid,ids);
334  load_crs_v<y,1,y,tgrad_u>(gradU_py,grid,ids);
335  load_crs_v<z,1,z,tgrad_u>(gradU_pz,grid,ids);
336 
337  load_crs_v<x,0,x,tgrad_u>(gradU_x,grid,ids);
338  load_crs_v<x,0,y,tgrad_u>(gradU_y,grid,ids);
339  load_crs_v<x,0,z,tgrad_u>(gradU_z,grid,ids);
340 
341  load_crs_v<x,1,x,tgrad_v>(gradV_px,grid,ids);
342  load_crs_v<y,1,y,tgrad_v>(gradV_py,grid,ids);
343  load_crs_v<z,1,z,tgrad_v>(gradV_pz,grid,ids);
344 
345  load_crs_v<x,0,x,tgrad_v>(gradV_x,grid,ids);
346  load_crs_v<x,0,y,tgrad_v>(gradV_y,grid,ids);
347  load_crs_v<x,0,z,tgrad_v>(gradV_z,grid,ids);
348 
349  load_crs<x,0,U>(Uc,grid,ids);
350  load_crs<x,0,V>(Vc,grid,ids);
351 
352  lapU += (gradU_px - gradU_x) / spacing[0];
353  lapV += (gradV_px - gradV_x) / spacing[0];
354  lapU += (gradU_py - gradU_y) / spacing[1];
355  lapV += (gradV_py - gradV_y) / spacing[1];
356  lapU += (gradU_pz - gradU_z) / spacing[2];
357  lapV += (gradV_pz - gradV_z) / spacing[2];
358 
359  // update based on Eq 2
360  outU = Uc + uFactor * lapU +
361  - deltaT * Uc * Vc * Vc +
362  - deltaT * F * (Uc - 1.0);
363 
364 
365  // update based on Eq 2
366  outV = Vc + vFactor * lapV +
367  deltaT * Uc * Vc * Vc +
368  - deltaT * (F+K) * Vc;
369 
370  Vc::Mask<double> surround;
371 
372  for (int i = 0 ; i < Vc::double_v::Size ; i++)
373  {surround[i] = (mask_sum[i] == 6);}
374 
375 
376  outU = Vc::iif(surround,outU,Uc);
377  outV = Vc::iif(surround,outV,Vc);
378 
379  store_crs<U_next>(grid,outU,ids);
380  store_crs<V_next>(grid,outV,ids);
381  };
382 
383  grid.template conv_cross_ids<1,double>({0,0,0},{sz[0]-1,sz[1] - 1,sz[2] - 1},func_lap);
384 
385 // New.write_frame("update",i);
386 
387  // Extend
388 
389  if (i % 5 == 0)
390  {
391  for (int j = 0 ; j < 2 ; j++)
392  {
393  if (j % 2 == 0)
394  {extend<U_next,V_next,U,V>(grid,sz,spacing);}
395  else
396  {extend<U,V,U_next,V_next>(grid,sz,spacing);}
397 
398  // Here we copy New into the old grid in preparation of the new step
399  // It would be better to alternate, but using this we can show the usage
400  // of the function copy. To note that copy work only on two grid of the same
401  // decomposition. If you want to copy also the decomposition, or force to be
402  // exactly the same, use Old = New
403  //New.copy_sparse(Old);
404  }
405  }
406 
407 /* auto it = grid.getDomainIterator();
408 
409  while (it.isNext())
410  {
411  // center point
412  auto Cp = it.get();
413 
414  // update based on Eq 2
415  grid.insert<U>(Cp) = grid.get<U_next>(Cp);
416  grid.insert<V>(Cp) = grid.get<V_next>(Cp);
417 
418  ++it;
419  }*/
420 
422 
423  // After copy we synchronize again the ghost part U and V
424  grid.ghost_get<U,V>();
425 
426  // Every 500 time step we output the configuration for
427  // visualization
428  if (i % 500 == 0)
429  {
430 // grid.save("output_" + std::to_string(count));
431  count++;
432  }
433 
434  if (v_cl.rank() == 0)
435  {std::cout << "STEP: " << i << " " << std::endl;}
436  if (i % 1000 == 0)
437  {
438  grid.write_frame("out",i);
439  }
440  }
441 
442  tot_sim.stop();
443  std::cout << "Total simulation: " << tot_sim.getwct() << std::endl;
444 
446 
458 
460  openfpm_finalize();
461 
463 
472 }
Point
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:27
timer::getwct
double getwct()
Return the elapsed real time.
Definition: timer.hpp:130
Ghost
Definition: Ghost.hpp:39
Point::norm
__device__ __host__ T norm()
norm of the vector
Definition: Point.hpp:231
grid_dist_id
This is a distributed grid.
Definition: grid_dist_id.hpp:96
Point::get
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition: Point.hpp:172
timer::start
void start()
Start the timer.
Definition: timer.hpp:90
Box
This class represent an N-dimensional box.
Definition: Box.hpp:60
Sphere
This class implement the Sphere concept in an N-dimensional space.
Definition: Sphere.hpp:23
cub::init
OutputIteratorT OffsetT ReductionOpT OuputT init
< [in] The initial value of the reduction
Definition: dispatch_reduce.cuh:119
periodicity
Boundary conditions.
Definition: common.hpp:21
timer
Class for cpu time benchmarking.
Definition: timer.hpp:27
timer::stop
void stop()
Stop the timer.
Definition: timer.hpp:119
F
[v_transform metafunction]
Definition: variadic_to_vmpl_unit_test.hpp:17