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
5#define FORTRAN_UPDATE
6
7#ifndef FORTRAN_UPDATE
8#include "Vc/Vc"
9#endif
10
48
49
50constexpr int x = 0;
51constexpr int y = 1;
52constexpr int z = 2;
53
54extern "C" void update_new(const int* lo, const int* hi,
55 double* u, const int* ulo, const int* uhi,
56 double* v, const int* vlo, const int* vhi,
57 double* flu, const int* fulo, const int* fuhi,
58 double* flv, const int* fvlo, const int* fvhi,
59 const double * dt, const double * uFactor, const double * vFactor, const double * F,
60 const double * K);
61
62
64
65void init(grid_dist_id<3,double,aggregate<double> > & OldU,
66 grid_dist_id<3,double,aggregate<double> > & OldV,
67 grid_dist_id<3,double,aggregate<double> > & NewU,
68 grid_dist_id<3,double,aggregate<double> > & NewV,
69 Box<3,double> & domain)
70{
71 auto it = OldU.getDomainIterator();
72
73 while (it.isNext())
74 {
75 // Get the local grid key
76 auto key = it.get();
77
78 // Old values U and V
79 OldU.get<0>(key) = 1.0;
80 OldV.get<0>(key) = 0.0;
81
82 // Old values U and V
83 NewU.get<0>(key) = 0.0;
84 NewV.get<0>(key) = 0.0;
85
86 ++it;
87 }
88
89 long int x_start = OldU.size(0)*1.55f/domain.getHigh(0);
90 long int y_start = OldU.size(1)*1.55f/domain.getHigh(1);
91 long int z_start = OldU.size(1)*1.55f/domain.getHigh(2);
92
93 long int x_stop = OldU.size(0)*1.85f/domain.getHigh(0);
94 long int y_stop = OldU.size(1)*1.85f/domain.getHigh(1);
95 long int z_stop = OldU.size(1)*1.85f/domain.getHigh(2);
96
97 grid_key_dx<3> start({x_start,y_start,z_start});
98 grid_key_dx<3> stop ({x_stop,y_stop,z_stop});
99 auto it_init = OldU.getSubDomainIterator(start,stop);
100
101 while (it_init.isNext())
102 {
103 auto key = it_init.get();
104
105 OldU.get<0>(key) = 0.5 + (((double)std::rand())/RAND_MAX -0.5)/10.0;
106 OldV.get<0>(key) = 0.25 + (((double)std::rand())/RAND_MAX -0.5)/20.0;
107
108 ++it_init;
109 }
110}
111
112
114
115void step(grid_dist_id<3, double, aggregate<double>> & OldU,
116 grid_dist_id<3, double, aggregate<double>> & OldV,
117 grid_dist_id<3, double, aggregate<double>> & NewU,
118 grid_dist_id<3, double, aggregate<double>> & NewV,
119 grid_key_dx<3> (& star_stencil_3D)[7],
120 double uFactor_s, double vFactor_s, double deltaT, double F, double K)
121{
122
123#ifndef FORTRAN_UPDATE
124
126
127 Vc::double_v uFactor = uFactor_s;
128 Vc::double_v vFactor = vFactor_s;
129
130 WHILE_M(OldU,star_stencil_3D)
131 auto & U_old = GET_GRID_M(OldU);
132 auto & V_old = GET_GRID_M(OldV);
133
134 auto & U_new = GET_GRID_M(NewU);
135 auto & V_new = GET_GRID_M(NewV);
136 ITERATE_3D_M(Vc::double_v::Size)
137
138 // center point
139 auto Cp = it.getStencil<0>();
140
141 // plus,minus X,Y,Z
142 auto mx = it.getStencil<1>();
143 auto px = it.getStencil<2>();
144 auto my = it.getStencil<3>();
145 auto py = it.getStencil<4>();
146 auto mz = it.getStencil<5>();
147 auto pz = it.getStencil<6>();
148
149 //
150 Vc::double_v u_c(&U_old.get<0>(Cp),Vc::Unaligned);
151 Vc::double_v u_mz(&U_old.get<0>(mz),Vc::Unaligned);
152 Vc::double_v u_pz(&U_old.get<0>(pz),Vc::Unaligned);
153 Vc::double_v u_my(&U_old.get<0>(my),Vc::Unaligned);
154 Vc::double_v u_py(&U_old.get<0>(py),Vc::Unaligned);
155 Vc::double_v u_mx(&U_old.get<0>(mx),Vc::Unaligned);
156 Vc::double_v u_px(&U_old.get<0>(px),Vc::Unaligned);
157
158
159 Vc::double_v v_c(&V_old.get<0>(Cp),Vc::Unaligned);
160 Vc::double_v v_mz(&V_old.get<0>(mz),Vc::Unaligned);
161 Vc::double_v v_pz(&V_old.get<0>(pz),Vc::Unaligned);
162 Vc::double_v v_my(&V_old.get<0>(my),Vc::Unaligned);
163 Vc::double_v v_py(&V_old.get<0>(py),Vc::Unaligned);
164 Vc::double_v v_mx(&V_old.get<0>(mx),Vc::Unaligned);
165 Vc::double_v v_px(&V_old.get<0>(px),Vc::Unaligned);
166
167 Vc::double_v out1 = u_c + uFactor * (u_mz + u_pz +
168 u_my + u_py +
169 u_mx + u_px +
170 - 6.0 * u_c) +
171 - deltaT * u_c * v_c * v_c
172 - deltaT * F * (u_c - 1.0);
173
174 Vc::double_v out2 = v_c + vFactor * (v_mz + v_pz +
175 v_my + v_py +
176 v_mx + v_px +
177 - 6.0 * v_c ) +
178 deltaT * u_c * v_c * v_c +
179 - deltaT * (F+K) * v_c;
180
181 out1.store(&U_new.get<0>(Cp),Vc::Unaligned);
182 out2.store(&V_new.get<0>(Cp),Vc::Unaligned);
183 END_LOOP_M(Vc::double_v::Size)
184
185
186
187#else
188
190
191 auto & ginfo = OldU.getLocalGridsInfo();
192
193 for (size_t i = 0 ; i < OldU.getN_loc_grid() ; i++)
194 {
195 auto & U_old = OldU.get_loc_grid(i);
196 auto & V_old = OldV.get_loc_grid(i);
197
198 auto & U_new = NewU.get_loc_grid(i);
199 auto & V_new = NewV.get_loc_grid(i);
200
201 int lo[3] = {(int)ginfo.get(i).Dbox.getLow(0),(int)ginfo.get(i).Dbox.getLow(1),(int)ginfo.get(i).Dbox.getLow(2)};
202 int hi[3] = {(int)ginfo.get(i).Dbox.getHigh(0),(int)ginfo.get(i).Dbox.getHigh(1),(int)ginfo.get(i).Dbox.getHigh(2)};
203
204 int ulo[3] = {0,0,0};
205 int uhi[3] = {(int)ginfo.get(i).GDbox.getHigh(0),(int)ginfo.get(i).GDbox.getHigh(1),(int)ginfo.get(i).GDbox.getHigh(2)};
206 int nulo[3] = {0,0,0};
207 int nuhi[3] = {(int)ginfo.get(i).GDbox.getHigh(0),(int)ginfo.get(i).GDbox.getHigh(1),(int)ginfo.get(i).GDbox.getHigh(2)};
208
209 int vlo[3] = {0,0,0};
210 int vhi[3] = {(int)ginfo.get(i).GDbox.getHigh(0),(int)ginfo.get(i).GDbox.getHigh(1),(int)ginfo.get(i).GDbox.getHigh(2)};
211 int nvlo[3] = {0,0,0};
212 int nvhi[3] = {(int)ginfo.get(i).GDbox.getHigh(0),(int)ginfo.get(i).GDbox.getHigh(1),(int)ginfo.get(i).GDbox.getHigh(2)};
213
214 update_new(lo,hi,
215 (double *)U_old.getPointer(),ulo,uhi,
216 (double *)V_old.getPointer(),vlo,vhi,
217 (double *)U_new.getPointer(),nulo,nuhi,
218 (double *)V_new.getPointer(),nulo,nvhi,
219 &deltaT, &uFactor_s, &vFactor_s,&F,&K);
220 }
221
223
224#endif
225}
226
228
229int main(int argc, char* argv[])
230{
231 openfpm_init(&argc,&argv);
232
233 // domain
234 Box<3,double> domain({0.0,0.0},{2.5,2.5,2.5});
235
236 // grid size
237 size_t sz[3] = {128,128,128};
238
239 // Define periodicity of the grid
240 periodicity<3> bc = {PERIODIC,PERIODIC,PERIODIC};
241
242 // Ghost in grid unit
244
245 // deltaT
246 double deltaT = 1.0;
247
248 // Diffusion constant for specie U
249 double du = 2*1e-5;
250
251 // Diffusion constant for specie V
252 double dv = 1*1e-5;
253
254 // Number of timesteps
255 size_t timeSteps = 5000;
256
257 // K and F (Physical constant in the equation)
258 double K = 0.053;
259 double F = 0.014;
260
262
275
276 grid_dist_id<3, double, aggregate<double>> OldU(sz,domain,g,bc);
277 grid_dist_id<3, double, aggregate<double>> OldV(OldU.getDecomposition(),sz,g);
278
279 // New grid with the decomposition of the old grid
280 grid_dist_id<3, double, aggregate<double>> NewU(OldU.getDecomposition(),sz,g);
281 grid_dist_id<3, double, aggregate<double>> NewV(OldV.getDecomposition(),sz,g);
282
283 // spacing of the grid on x and y
284
285 double spacing[3] = {OldU.spacing(0),OldU.spacing(1),OldU.spacing(2)};
286
287 init(OldU,OldV,NewU,NewV,domain);
288
290
291 // sync the ghost
292 size_t count = 0;
293 OldU.template ghost_get<0>();
294 OldV.template ghost_get<0>();
295
296 // because we assume that spacing[x] == spacing[y] we use formula 2
297 // and we calculate the prefactor of Eq 2
298 double uFactor = deltaT * du/(spacing[x]*spacing[x]);
299 double vFactor = deltaT * dv/(spacing[x]*spacing[x]);
300
301 timer tot_sim;
302 tot_sim.start();
303
304 auto & v_cl = create_vcluster();
305
306 static grid_key_dx<3> star_stencil_3D[7] = {{0,0,0},
307 {0,0,-1},
308 {0,0,1},
309 {0,-1,0},
310 {0,1,0},
311 {-1,0,0},
312 {1,0,0}};
313
314 for (size_t i = 0; i < timeSteps; ++i)
315 {
316 if (i % 100 == 0 && v_cl.rank() == 0)
317 {
318 std::cout << "STEP: " << i << std::endl;
319 }
320
353
354 if (i % 2 == 0)
355 {
356 step(OldU,OldV,NewU,NewV,star_stencil_3D,uFactor,vFactor,deltaT,F,K);
357
358 NewU.ghost_get<0>();
359 NewV.ghost_get<0>();
360 }
361 else
362 {
363 step(NewU,NewV,OldU,OldV,star_stencil_3D,uFactor,vFactor,deltaT,F,K);
364
365 OldU.ghost_get<0>();
366 OldV.ghost_get<0>();
367 }
368
370
381
382 // Every 2000 time step we output the configuration on hdf5
383/* if (i % 2000 == 0)
384 {
385 OldU.save("output_u_" + std::to_string(count));
386 OldV.save("output_v_" + std::to_string(count));
387 count++;
388 }*/
389
391 }
392
393 tot_sim.stop();
394
395 if (create_vcluster().rank() == 0)
396 {std::cout << "Total simulation: " << tot_sim.getwct() << std::endl;}
397
398 // We frite the final configuration
399 OldV.write("final");
400
402
415
416 openfpm_finalize();
417
419
428}
429
430
This class represent an N-dimensional box.
Definition Box.hpp:61
This is a distributed grid.
grid_key_dx is the key to access any element in the grid
Definition grid_key.hpp:19
__device__ __host__ index_type get(index_type i) const
Get the i index.
Definition grid_key.hpp:503
Class for cpu time benchmarking.
Definition timer.hpp:28
void stop()
Stop the timer.
Definition timer.hpp:119
void start()
Start the timer.
Definition timer.hpp:90
double getwct()
Return the elapsed real time.
Definition timer.hpp:130
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]
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Boundary conditions.
Definition common.hpp:22