OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
main.cu
1//#define VCLUSTER_PERF_REPORT
2//#define SYNC_BEFORE_TAKE_TIME
3//#define ENABLE_GRID_DIST_ID_PERF_STATS
4#include "Grid/grid_dist_id.hpp"
5#include "data_type/aggregate.hpp"
6#include "timer.hpp"
7
38#ifdef __NVCC__
39
40constexpr int U = 0;
41constexpr int V = 1;
42constexpr int U_next = 2;
43constexpr int V_next = 3;
44
45typedef sgrid_dist_id_gpu<3,double,aggregate<double,double,double,double> > sgrid_type;
46
47void init(sgrid_type & grid, Box<3,double> & domain)
48{
49 auto it = grid.getGridIterator();
50 Point<3,double> p[8]= {{0.35,0.35,0.35},
51 {0.35,2.0,2.0},
52 {2.0,0.35,2.0},
53 {2.0,2.0,0.35},
54 {0.35,0.35,2.0},
55 {0.35,2.0,0.35},
56 {2.0,0.35,0.35},
57 {2.0,2.0,2.0}};
58
59
60// Point<3,double> u({1.0,0.0,0.0});
61// Box<3,double> channel_box(p3,p1);
62
63 double spacing_x = grid.spacing(0);
64 double spacing_y = grid.spacing(1);
65 double spacing_z = grid.spacing(2);
66
67 typedef typename GetAddBlockType<sgrid_type>::type InsertBlockT;
68
69 // Draw spheres
70 for (int i = 0 ; i < 8 ; i++)
71 {
72 Sphere<3,double> sph(p[i],0.3);
73
74 Box<3,size_t> bx;
75
76 for (int i = 0 ; i < 3 ; i++)
77 {
78 bx.setLow(i,(size_t)((sph.center(i) - 0.31)/grid.spacing(i)));
79 bx.setHigh(i,(size_t)((sph.center(i) + 0.31)/grid.spacing(i)));
80 }
81
82 grid.addPoints(bx.getKP1(),bx.getKP2(),[spacing_x,spacing_y,spacing_z,sph] __device__ (int i, int j, int k)
83 {
84 Point<3,double> pc({i*spacing_x,j*spacing_y,k*spacing_z});
85
86 // Check if the point is in the domain
87 if (sph.isInside(pc) )
88 {return true;}
89
90 return false;
91 },
92 [] __device__ (InsertBlockT & data, int i, int j, int k)
93 {
94 data.template get<U>() = 1.0;
95 data.template get<V>() = 0.0;
96 }
97 );
98
99 grid.template flush<smax_<U>,smax_<V>>(flush_type::FLUSH_ON_DEVICE);
100 grid.removeUnusedBuffers();
101 }
102
103 //channels
104
105 Box<3,double> b({0.25,0.25,0.25},{2.1,2.1,2.1});
106
107 for (int k = 0 ; k < 3 ; k++)
108 {
109 for (int s = 0 ; s < 2 ; s++)
110 {
111 for (int i = 0 ; i < 2 ; i++)
112 {
113 Point<3,double> u({1.0*(((s+i)%2) == 0 && k != 2),1.0*(((s+i+1)%2) == 0 && k != 2),(k == 2)*1.0});
114 Point<3,double> c({(i == 0)?0.35:2.0,(s == 0)?0.35:2.0,(k == 0)?0.35:2.0});
115
116 Box<3,size_t> bx;
117
118 for (int i = 0 ; i < 3 ; i++)
119 {
120 if (c[i] == 2.0)
121 {
122 if (u[i] == 1.0)
123 {
124 bx.setLow(i,(size_t)(0.34/grid.spacing(i)));
125 bx.setHigh(i,(size_t)(2.01/grid.spacing(i)));
126 }
127 else
128 {
129 bx.setLow(i,(size_t)((c[i] - 0.11)/grid.spacing(i)));
130 bx.setHigh(i,(size_t)((c[i] + 0.11)/grid.spacing(i)));
131 }
132 }
133 else
134 {
135 if (u[i] == 1.0)
136 {
137 bx.setLow(i,(size_t)(0.34/grid.spacing(i)));
138 bx.setHigh(i,(size_t)(2.01/grid.spacing(i)));
139 }
140 else
141 {
142 bx.setLow(i,(size_t)((c[i] - 0.11)/grid.spacing(i)));
143 bx.setHigh(i,(size_t)((c[i] + 0.11)/grid.spacing(i)));
144 }
145 }
146 }
147
148 grid.addPoints(bx.getKP1(),bx.getKP2(),[spacing_x,spacing_y,spacing_z,u,c,b] __device__ (int i, int j, int k)
149 {
150 Point<3,double> pc({i*spacing_x,j*spacing_y,k*spacing_z});
151 Point<3,double> pcs({i*spacing_x,j*spacing_y,k*spacing_z});
152 Point<3,double> vp;
153
154 // shift
155 pc -= c;
156
157 // calculate the distance from the diagonal
158 vp.get(0) = pc.get(1)*u.get(2) - pc.get(2)*u.get(1);
159 vp.get(1) = pc.get(2)*u.get(0) - pc.get(0)*u.get(2);
160 vp.get(2) = pc.get(0)*u.get(1) - pc.get(1)*u.get(0);
161
162 double distance = vp.norm();
163
164 // Check if the point is in the domain
165 if (distance < 0.1 && b.isInside(pcs) == true )
166 {return true;}
167
168 return false;
169 },
170 [] __device__ (InsertBlockT & data, int i, int j, int k)
171 {
172 data.template get<U>() = 1.0;
173 data.template get<V>() = 0.0;
174 }
175 );
176
177 grid.template flush<smax_<U>,smax_<V>>(flush_type::FLUSH_ON_DEVICE);
178 grid.removeUnusedBuffers();
179 }
180 }
181 }
182
183 // cross channel
184
185 int s = 0;
186 for (int s = 0 ; s < 2 ; s++)
187 {
188 for (int i = 0 ; i < 2 ; i++)
189 {
190 Point<3,double> c({(i == 0)?0.35:2.0,(s == 0)?0.35:2.0,0.35});
191 Point<3,double> u({(i == 0)?1.0:-1.0,(s == 0)?1.0:-1.0,1.0});
192
193 Box<3,size_t> bx;
194
195 for (int k = 0 ; k < 16; k++)
196 {
197 for (int s = 0 ; s < 3 ; s++)
198 {
199 if (u[s] > 0.0)
200 {
201 bx.setLow(s,(c[s] + k*(u[s]/9.0))/grid.spacing(s) );
202 bx.setHigh(s,(c[s] + (k+3)*(u[s]/9.0) )/ grid.spacing(s) );
203 }
204 else
205 {
206 bx.setLow(s,(c[s] + (k+3)*(u[s]/9.0) )/grid.spacing(s) );
207 bx.setHigh(s,(c[s] + k*(u[s]/9.0))/ grid.spacing(s) );
208 }
209 }
210
211 grid.addPoints(bx.getKP1(),bx.getKP2(),[spacing_x,spacing_y,spacing_z,u,c,b] __device__ (int i, int j, int k)
212 {
213 Point<3,double> pc({i*spacing_x,j*spacing_y,k*spacing_z});
214 Point<3,double> pcs({i*spacing_x,j*spacing_y,k*spacing_z});
215 Point<3,double> vp;
216
217 // shift
218 pc -= c;
219
220 // calculate the distance from the diagonal
221 vp.get(0) = pc.get(1)*u.get(2) - pc.get(2)*u.get(1);
222 vp.get(1) = pc.get(2)*u.get(0) - pc.get(0)*u.get(2);
223 vp.get(2) = pc.get(0)*u.get(1) - pc.get(1)*u.get(0);
224
225 double distance = vp.norm() / sqrt(3.0);
226
227 // Check if the point is in the domain
228 if (distance < 0.1 && b.isInside(pcs) == true )
229 {return true;}
230
231 return false;
232 },
233 [] __device__ (InsertBlockT & data, int i, int j, int k)
234 {
235 data.template get<U>() = 1.0;
236 data.template get<V>() = 0.0;
237 }
238 );
239
240 grid.template flush<smax_<U>,smax_<V>>(flush_type::FLUSH_ON_DEVICE);
241 grid.removeUnusedBuffers();
242 }
243 }
244
245 }
246
247 long int x_start = grid.size(0)*1.95f/domain.getHigh(0);
248 long int y_start = grid.size(1)*1.95f/domain.getHigh(1);
249 long int z_start = grid.size(1)*1.95f/domain.getHigh(2);
250
251 long int x_stop = grid.size(0)*2.05f/domain.getHigh(0);
252 long int y_stop = grid.size(1)*2.05f/domain.getHigh(1);
253 long int z_stop = grid.size(1)*2.05f/domain.getHigh(2);
254
255 grid_key_dx<3> start({x_start,y_start,z_start});
256 grid_key_dx<3> stop ({x_stop,y_stop,z_stop});
257
258 grid.addPoints(start,stop,[] __device__ (int i, int j, int k)
259 {
260 return true;
261 },
262 [] __device__ (InsertBlockT & data, int i, int j, int k)
263 {
264 data.template get<U>() = 0.5;
265 data.template get<V>() = 0.24;
266 }
267 );
268
269 grid.template flush<smin_<U>,smax_<V>>(flush_type::FLUSH_ON_DEVICE);
270
271 grid.removeUnusedBuffers();
272}
273
274
275int main(int argc, char* argv[])
276{
277 openfpm_init(&argc,&argv);
278
279 // domain
280 Box<3,double> domain({0.0,0.0,0.0},{2.5,2.5,2.5});
281
282 // grid size
283 size_t sz[3] = {384,384,384};
284
285 // Define periodicity of the grid
286 periodicity<3> bc = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
287
288 // Ghost in grid unit
289 Ghost<3,long int> g(1);
290
291 // deltaT
292 double deltaT = 0.2;
293
294 // Diffusion constant for specie U
295 double du = 2*1e-5;
296
297 // Diffusion constant for specie V
298 double dv = 1*1e-5;
299
300#ifdef TEST_RUN
301 // Number of timesteps
302 size_t timeSteps = 300;
303#else
304 // Number of timesteps
305 size_t timeSteps = 50000;
306#endif
307
308 // K and F (Physical constant in the equation)
309 double K = 0.053;
310 double F = 0.014;
311
312 sgrid_type grid(sz,domain,g,bc);
313
314 grid.template setBackgroundValue<0>(-0.5);
315 grid.template setBackgroundValue<1>(-0.5);
316 grid.template setBackgroundValue<2>(-0.5);
317 grid.template setBackgroundValue<3>(-0.5);
318
319 // spacing of the grid on x and y
320 double spacing[3] = {grid.spacing(0),grid.spacing(1),grid.spacing(2)};
321
322 init(grid,domain);
323
324 // sync the ghost
325 grid.template ghost_get<U,V>(RUN_ON_DEVICE);
326
327 // because we assume that spacing[x] == spacing[y] we use formula 2
328 // and we calculate the prefactor of Eq 2
329 double uFactor = deltaT * du/(spacing[0]*spacing[0]);
330 double vFactor = deltaT * dv/(spacing[0]*spacing[0]);
331
332 grid.template deviceToHost<U,V>();
333
334 timer tot_sim;
335 tot_sim.start();
336
337 for (size_t i = 0; i < timeSteps; ++i)
338 {
340
341 typedef typename GetCpBlockType<decltype(grid),0,1>::type CpBlockType;
342
343 auto func = [uFactor,vFactor,deltaT,F,K] __device__ (double & u_out, double & v_out,
344 CpBlockType & u, CpBlockType & v,
345 int i, int j, int k){
346
347 double uc = u(i,j,k);
348 double vc = v(i,j,k);
349
350 double u_px = u(i+1,j,k);
351 double u_mx = u(i-1,j,k);
352
353 double u_py = u(i,j+1,k);
354 double u_my = u(i,j-1,k);
355
356 double u_pz = u(i,j,k+1);
357 double u_mz = u(i,j,k-1);
358
359 double v_px = v(i+1,j,k);
360 double v_mx = v(i-1,j,k);
361
362 double v_py = v(i,j+1,k);
363 double v_my = v(i,j-1,k);
364
365 double v_pz = v(i,j,k+1);
366 double v_mz = v(i,j,k-1);
367
368 // U fix
369
370 if (u_mx < -0.1 && u_px < -0.1)
371 {
372 u_mx = uc;
373 u_px = uc;
374 }
375
376 if (u_mx < -0.1)
377 {u_mx = u_px;}
378
379 if (u_px < -0.1)
380 {u_px = u_mx;}
381
382 if (u_my < -0.1 && u_py < -0.1)
383 {
384 u_my = uc;
385 u_py = uc;
386 }
387
388 if (u_my < -0.1)
389 {u_my = u_py;}
390
391 if (u_py < -0.1)
392 {u_py = u_my;}
393
394 if (u_mz < -0.1 && u_pz < -0.1)
395 {
396 u_mz = uc;
397 u_pz = uc;
398 }
399
400 if (u_mz < -0.1)
401 {u_mz = u_pz;}
402
403 if (u_pz < -0.1)
404 {u_pz = u_mz;}
405
406 // V fix
407
408 if (v_mx < -0.1 && v_px < -0.1)
409 {
410 v_mx = uc;
411 v_px = uc;
412 }
413
414 if (v_mx < -0.1)
415 {v_mx = v_px;}
416
417 if (v_px < -0.1)
418 {v_px = v_mx;}
419
420 if (v_my < -0.1 && v_py < -0.1)
421 {
422 v_my = uc;
423 v_py = uc;
424 }
425
426 if (v_my < -0.1)
427 {v_my = v_py;}
428
429 if (v_py < -0.1)
430 {v_py = v_my;}
431
432 if (v_mz < -0.1 && v_pz < -0.1)
433 {
434 v_mz = uc;
435 v_pz = uc;
436 }
437
438 if (v_mz < -0.1)
439 {v_mz = v_pz;}
440
441 if (v_pz < -0.1)
442 {v_pz = v_mz;}
443
444 u_out = uc + uFactor *(u_mx + u_px +
445 u_my + u_py +
446 u_mz + u_pz - 6.0*uc) - deltaT * uc*vc*vc
447 - deltaT * F * (uc - 1.0);
448
449
450 v_out = vc + vFactor *(v_mx + v_px +
451 v_py + v_my +
452 v_mz + v_pz - 6.0*vc) + deltaT * uc*vc*vc
453 - deltaT * (F+K) * vc;
454
455 };
456
457 if (i % 2 == 0)
458 {
459 grid.conv2<U,V,U_next,V_next,1>({0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func);
460
461 cudaDeviceSynchronize();
462
463 // After copy we synchronize again the ghost part U and V
464
465 grid.ghost_get<U_next,V_next>(RUN_ON_DEVICE | SKIP_LABELLING);
466 }
467 else
468 {
469 grid.conv2<U_next,V_next,U,V,1>({0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func);
470
471 cudaDeviceSynchronize();
472
473 // After copy we synchronize again the ghost part U and V
474 grid.ghost_get<U,V>(RUN_ON_DEVICE | SKIP_LABELLING);
475 }
476
478
479 // After copy we synchronize again the ghost part U and V
480
481 // Every 500 time step we output the configuration for
482 // visualization
483/* if (i % 500 == 0)
484 {
485 grid.save("output_" + std::to_string(count));
486 count++;
487 }*/
488
489 std::cout << "STEP: " << i << std::endl;
490/* if (i % 300 == 0)
491 {
492 grid.template deviceToHost<U,V>();
493 grid.write_frame("out",i);
494 }*/
495 }
496
497 tot_sim.stop();
498 std::cout << "Total simulation: " << tot_sim.getwct() << std::endl;
499
500 grid.print_stats();
501
502 create_vcluster().print_stats();
503
504 grid.template deviceToHost<U,V>();
505 grid.write("Final");
506
508
521
522 openfpm_finalize();
523
525
534}
535
536#else
537
538int main(int argc, char* argv[])
539{
540 return 0;
541}
542
543#endif
544
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
Box::setHigh
__device__ __host__ void setHigh(int i, T val)
set the high interval of the box
Definition Box.hpp:544
Box::getKP2
grid_key_dx< dim > getKP2() const
Get the point p12 as grid_key_dx.
Definition Box.hpp:669
Box::setLow
__device__ __host__ void setLow(int i, T val)
set the low interval of the box
Definition Box.hpp:533
Box::getKP1
grid_key_dx< dim > getKP1() const
Get the point p1 as grid_key_dx.
Definition Box.hpp:656
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
grid_key_dx
grid_key_dx is the key to access any element in the grid
Definition grid_key.hpp:19
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::int
KeyT const ValueT ValueT OffsetIteratorT OffsetIteratorT int
[in] The number of segments that comprise the sorting data
Definition dispatch_radix_sort.cuh:336
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
GetCpBlockType
get the type of the block
Definition SparseGridGpu.hpp:162
enc_num
Definition encap_num.hpp:13
periodicity
Boundary conditions.
Definition common.hpp:22
smax_
Definition map_vector_sparse_cuda_kernels.cuh:192