OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
main.cpp
1#include "Grid/grid_dist_id.hpp"
2#include "data_type/aggregate.hpp"
3#include "timer.hpp"
4
9constexpr int U = 0;
10constexpr int V = 1;
11constexpr int phi = 2;
12constexpr int normal = 3;
13constexpr int tgrad_u = 4;
14constexpr int tgrad_v = 5;
15constexpr int U_next = 6;
16constexpr int V_next = 7;
17
18constexpr int x = 0;
19constexpr int y = 1;
20constexpr int z = 2;
21
22typedef sgrid_dist_soa<3,double,aggregate<double,double,double,double[3],double[3],double[3],double,double> > sgrid_type;
23
24void 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
87template<unsigned int U_src,unsigned int V_src,unsigned int U_dst, unsigned int V_dst>
88void 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
180int 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
459
460 openfpm_finalize();
461
463
472}
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
Point
This class implement the point shape in an N-dimensional space.
Definition Point.hpp:28
Point::get
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition Point.hpp:172
Point::norm
__device__ __host__ T norm() const
norm of the vector
Definition Point.hpp:231
Sphere
This class implement the Sphere concept in an N-dimensional space.
Definition Sphere.hpp:24
grid_dist_id
This is a distributed grid.
Definition grid_dist_id.hpp:278
timer
Class for cpu time benchmarking.
Definition timer.hpp:28
timer::stop
void stop()
Stop the timer.
Definition timer.hpp:119
timer::start
void start()
Start the timer.
Definition timer.hpp:90
timer::getwct
double getwct()
Return the elapsed real time.
Definition timer.hpp:130
cub::init
OutputIteratorT OffsetT ReductionOpT OuputT init
< [in] The initial value of the reduction
Definition dispatch_reduce.cuh:119
F
[v_transform metafunction]
Definition variadic_to_vmpl_unit_test.hpp:18
periodicity
Boundary conditions.
Definition common.hpp:22