4#include "Grid/grid_dist_id.hpp"
5#include "data_type/aggregate.hpp"
42constexpr int U_next = 2;
43constexpr int V_next = 3;
45typedef sgrid_dist_id_gpu<3,double,aggregate<double,double,double,double> >
sgrid_type;
49 auto it =
grid.getGridIterator();
63 double spacing_x =
grid.spacing(0);
64 double spacing_y =
grid.spacing(1);
65 double spacing_z =
grid.spacing(2);
70 for (
int i = 0 ; i < 8 ; i++)
76 for (
int i = 0 ; i < 3 ; i++)
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)));
82 grid.addPoints(bx.
getKP1(),bx.
getKP2(),[spacing_x,spacing_y,spacing_z,sph] __device__ (
int i,
int j,
int k)
87 if (sph.isInside(pc) )
92 [] __device__ (InsertBlockT & data,
int i,
int j,
int k)
94 data.template get<U>() = 1.0;
95 data.template get<V>() = 0.0;
99 grid.template flush<smax_<U>,
smax_<V>>(flush_type::FLUSH_ON_DEVICE);
100 grid.removeUnusedBuffers();
107 for (
int k = 0 ; k < 3 ; k++)
109 for (
int s = 0 ; s < 2 ; s++)
111 for (
int i = 0 ; i < 2 ; i++)
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});
118 for (
int i = 0 ; i < 3 ; i++)
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)));
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)));
148 grid.addPoints(bx.
getKP1(),bx.
getKP2(),[spacing_x,spacing_y,spacing_z,u,c,b] __device__ (
int i,
int j,
int k)
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);
162 double distance = vp.
norm();
165 if (distance < 0.1 && b.isInside(pcs) ==
true )
170 [] __device__ (InsertBlockT & data,
int i,
int j,
int k)
172 data.template get<U>() = 1.0;
173 data.template get<V>() = 0.0;
177 grid.template flush<smax_<U>,
smax_<V>>(flush_type::FLUSH_ON_DEVICE);
178 grid.removeUnusedBuffers();
186 for (
int s = 0 ; s < 2 ; s++)
188 for (
int i = 0 ; i < 2 ; i++)
195 for (
int k = 0 ; k < 16; k++)
197 for (
int s = 0 ; s < 3 ; s++)
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) );
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) );
211 grid.addPoints(bx.
getKP1(),bx.
getKP2(),[spacing_x,spacing_y,spacing_z,u,c,b] __device__ (
int i,
int j,
int k)
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);
225 double distance = vp.
norm() / sqrt(3.0);
228 if (distance < 0.1 && b.isInside(pcs) ==
true )
233 [] __device__ (InsertBlockT & data,
int i,
int j,
int k)
235 data.template get<U>() = 1.0;
236 data.template get<V>() = 0.0;
240 grid.template flush<smax_<U>,
smax_<V>>(flush_type::FLUSH_ON_DEVICE);
241 grid.removeUnusedBuffers();
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);
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);
258 grid.addPoints(start,stop,[] __device__ (
int i,
int j,
int k)
262 [] __device__ (InsertBlockT & data,
int i,
int j,
int k)
264 data.template get<U>() = 0.5;
265 data.template get<V>() = 0.24;
269 grid.template flush<smin_<U>,
smax_<V>>(flush_type::FLUSH_ON_DEVICE);
271 grid.removeUnusedBuffers();
275int main(
int argc,
char* argv[])
277 openfpm_init(&argc,&argv);
283 size_t sz[3] = {384,384,384};
302 size_t timeSteps = 300;
305 size_t timeSteps = 50000;
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);
320 double spacing[3] = {
grid.spacing(0),
grid.spacing(1),
grid.spacing(2)};
325 grid.template ghost_get<U,V>(RUN_ON_DEVICE);
329 double uFactor = deltaT * du/(spacing[0]*spacing[0]);
330 double vFactor = deltaT * dv/(spacing[0]*spacing[0]);
332 grid.template deviceToHost<U,V>();
337 for (
size_t i = 0; i < timeSteps; ++i)
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){
347 double uc = u(i,j,k);
348 double vc = v(i,j,k);
350 double u_px = u(i+1,j,k);
351 double u_mx = u(i-1,j,k);
353 double u_py = u(i,j+1,k);
354 double u_my = u(i,j-1,k);
356 double u_pz = u(i,j,k+1);
357 double u_mz = u(i,j,k-1);
359 double v_px = v(i+1,j,k);
360 double v_mx = v(i-1,j,k);
362 double v_py = v(i,j+1,k);
363 double v_my = v(i,j-1,k);
365 double v_pz = v(i,j,k+1);
366 double v_mz = v(i,j,k-1);
370 if (u_mx < -0.1 && u_px < -0.1)
382 if (u_my < -0.1 && u_py < -0.1)
394 if (u_mz < -0.1 && u_pz < -0.1)
408 if (v_mx < -0.1 && v_px < -0.1)
420 if (v_my < -0.1 && v_py < -0.1)
432 if (v_mz < -0.1 && v_pz < -0.1)
444 u_out = uc + uFactor *(u_mx + u_px +
446 u_mz + u_pz - 6.0*uc) - deltaT * uc*vc*vc
447 - deltaT *
F * (uc - 1.0);
450 v_out = vc + vFactor *(v_mx + v_px +
452 v_mz + v_pz - 6.0*vc) + deltaT * uc*vc*vc
453 - deltaT * (
F+K) * vc;
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);
461 cudaDeviceSynchronize();
465 grid.ghost_get<U_next,V_next>(RUN_ON_DEVICE | SKIP_LABELLING);
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);
471 cudaDeviceSynchronize();
474 grid.ghost_get<U,V>(RUN_ON_DEVICE | SKIP_LABELLING);
489 std::cout <<
"STEP: " << i << std::endl;
498 std::cout <<
"Total simulation: " << tot_sim.
getwct() << std::endl;
502 create_vcluster().print_stats();
504 grid.template deviceToHost<U,V>();
538int main(
int argc,
char* argv[])
This class represent an N-dimensional box.
__device__ __host__ void setHigh(int i, T val)
set the high interval of the box
grid_key_dx< dim > getKP2() const
Get the point p12 as grid_key_dx.
__device__ __host__ void setLow(int i, T val)
set the low interval of the box
grid_key_dx< dim > getKP1() const
Get the point p1 as grid_key_dx.
This class implement the point shape in an N-dimensional space.
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
__device__ __host__ T norm() const
norm of the vector
This class implement the Sphere concept in an N-dimensional space.
This is a distributed grid.
grid_key_dx is the key to access any element in the grid
Class for cpu time benchmarking.
void stop()
Stop the timer.
void start()
Start the timer.
double getwct()
Return the elapsed real time.
KeyT const ValueT ValueT OffsetIteratorT OffsetIteratorT int
[in] The number of segments that comprise the sorting data
OutputIteratorT OffsetT ReductionOpT OuputT init
< [in] The initial value of the reduction
[v_transform metafunction]
get the type of the block